Front. Microbiol. Frontiers in Microbiology Front. Microbiol. 1664-302X Frontiers Media S.A. 10.3389/fmicb.2020.00129 Microbiology Original Research Broad and Fine Scale Variability in Bacterial Diversity and Cyanotoxin Quotas in Benthic Cyanobacterial Mats Thomson-Laing Georgia 1 * Puddick Jonathan 1 Laroche Olivier 1 2 Fulton Samantha 1 Steiner Konstanze 1 Heath Mark W. 3 Wood Susanna A. 1 1Cawthron Institute, Nelson, New Zealand 2Department of Oceanography, School of Ocean and Earth Science and Technology, University of Hawai‘i at Mānoa, Honolulu, HI, United States 3Greater Wellington Reginal Council, Wellington, New Zealand

Edited by: Jonathan P. Zehr, University of California, Santa Cruz, United States

Reviewed by: Nico Salmaso, Fondazione Edmund Mach, Italy; Anne Bernhard, Connecticut College, United States

*Correspondence: Georgia Thomson-Laing, Georgia.Thomson-Laing@cawthron.org.nz

This article was submitted to Aquatic Microbiology, a section of the journal Frontiers in Microbiology

06 02 2020 2020 11 129 27 05 2019 21 01 2020 Copyright © 2020 Thomson-Laing, Puddick, Laroche, Fulton, Steiner, Heath and Wood. 2020 Thomson-Laing, Puddick, Laroche, Fulton, Steiner, Heath and Wood

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Benthic proliferations of Microcoleus autumnalis (basionym Phormidium autumnale) and closely related taxa are being reported with increasing frequency in streams and rivers worldwide. This species commonly produces the potent neurotoxin anatoxin, and exposure to this has resulted in animal fatalities and human health concerns. Bacterial communities within cyanobacterial assemblages can facilitate processes such as nutrient cycling and are posited to influence cyanobacterial growth and function. However, there is limited knowledge on spatial variability of bacterial communities associated with benthic cyanobacteria and anatoxin content and quotas. In this study, M. autumnalis-dominated mat samples were collected from six sites in two New Zealand streams. Associated bacterial communities were characterized using 16S rRNA metabarcoding, anatoxin content by liquid chromatography-mass spectrometry and anaC copies using droplet digital PCR. Bacterial assemblages differed significantly when amplicon sequence variants were compared between streams and most sites within streams. These differences were associated with conductivity, DRP, DIN, temperature, anatoxin concentration, and quota. Despite the differences in bacterial community composition; at phyla, class and order levels there was high similarity across spatial scales, with Bacteroidetes (ca. 67%) and Proteobacteria (ca. 25%) dominant. There was significant variability in total anatoxin concentrations between sites in both streams (p < 0.001). When the data were converted to anatoxin quotas variability was reduced, suggesting that the relative abundance of toxic genotypes is a key driver of total anatoxin concentrations in mats. This study demonstrates the complexity of microbial communities within M. autumnalis-dominated mats and highlights their likely important role in within-mat nutrient cycling processes.

cyanobacteria Phormidium Microcoleus biofilm high-throughput sequencing 16S rRNA metabarcoding droplet digital PCR liquid chromatography-mass spectrometry

香京julia种子在线播放

    1. <form id=HxFbUHhlv><nobr id=HxFbUHhlv></nobr></form>
      <address id=HxFbUHhlv><nobr id=HxFbUHhlv><nobr id=HxFbUHhlv></nobr></nobr></address>

      Introduction

      An increasing number of cyanobacterial species are now known to produce cyanotoxins (Mowe et al., 2015; Cirés et al., 2017). Their proliferation, especially in freshwater environments, threatens both water quality and human and animal health (Codd et al., 1999; Codd, 2000; Paerl et al., 2001; Lee et al., 2017; Huisman et al., 2018). The intensification of anthropogenic impacts, as well as changing climates (O’Neil et al., 2012; Paerl and Paul, 2012; Huisman et al., 2018) has led to an increased occurrence of cyanobacterial blooms in both lentic and lotic environments. Research on harmful planktonic cyanobacterial blooms (primarily in lakes) has been undertaken for decades (e.g., Chorus and Bartram, 1999). Whilst toxic benthic cyanobacteria in streams and lakes pose similar issues, research in this area has only recently intensified (Quiblier et al., 2013).

      Recent research on toxic benthic cyanobacteria has shown the presence of many cyanotoxins in benthic cyanobacterial mats worldwide (reviewed by Quiblier et al., 2013). For example, microcystins have been identified in benthic samples from Spanish reservoirs (Aboal and Puig, 2005), alpine lakes in Switzerland (Mez et al., 1997), the Nile River and irrigation canals in Egypt (Mohamed et al., 2006), Brazilian streams and water supply reservoirs (Borges et al., 2015) and Californian streams (Fetscher et al., 2015). Other toxins including lyngbyatoxin, saxitoxins, cylindrospermopsin, and nodularin have also been identified in freshwater stream biofilms (Seifert et al., 2007; Borges et al., 2015; Fetscher et al., 2015; Foss et al., 2017, 2018). The most widely reported toxins from benthic mats are anatoxins, with production commonly attributed to Microcoleus or Phormidium. For example, anatoxins have been reported in mats from streams in France (in 2003; Gugger et al., 2005), New Zealand (in 2005; Wood et al., 2007), Netherlands (in 2011; Faassen et al., 2012), and in California, United States (2013–2015; Bouma-Gregson et al., 2018a). In these cases, research was undertaken in response to rapid dog deaths following consumption of toxic benthic cyanobacteria. Although there has been an increase in the identification of toxin-producing benthic cyanobacteria and their potential risks, most studies remain descriptive, i.e., describing toxin detection and causative species.

      In New Zealand rivers and streams the most prevalent benthic cyanobacteria is the filamentous Microcoleus autumnalis (basionym Phormidium autumnale), which commonly produces anatoxins. There has been an increase in reports of M. autumnalis-dominated proliferations (defined in New Zealand as greater than 20% mat coverage of stream bed substrate) in recent years (McAllister et al., 2016). Although dominated by M. autumnalis, these mats usually contain many other organisms (e.g., bacteria, diatoms and eukaryotic algae) and inorganic matter (e.g., sediment), bound together by extracellular polymeric substances (Hart et al., 2013; Brasell et al., 2015; Wood et al., 2015). Consortial interactions between microbial communities, specifically heterotrophic bacteria and cyanobacteria are believed to alleviate metabolic constraints through the exchange of carbon, specific growth factors (e.g., vitamins) and limiting nutrients (Paerl et al., 2000, 2001). Molecular approaches have now been extensively used for assessing the diversity of heterotrophic prokaryotes associated with cyanobacteria, with high throughput sequencing-based approaches now commonplace (Cai et al., 2014; Louati et al., 2015; Woodhouse et al., 2015; Berry et al., 2017; Parulekar et al., 2017; Su et al., 2017). To-date most studies have focused on planktonic cyanobacteria, with studies showing that cyanobacterial blooms disturb the natural bacterial community diversity and composition, but that the community returns relatively quickly after the blooms dissipate. Berry et al. (2017) found that Microcystis and Synechococcus bloom events in Lake Erie coincided with large shifts in non-cyanobacterial community composition that could also be predicted by environmental parameters such as pH, chlorophyll-a, temperature, and water movements. Some taxa (i.e., Actinobacteria) exhibited more pronounced changes and sensitivity to these shifts in bloom dynamics.

      Three recent studies have investigated relationships between M. autumnalis-dominated mats and associated microbial communities. Brasell et al. (2015) studied shifts in microbial communities through successional cycles of M. autumnalis-dominated mats and showed three distinct phases. Echenique-Subiabre et al. (2018) identified similar bacterial communities in Microcoleus-dominated mats from both New Zealand and French streams, with Bray-Curtis dissimilarities between French and New Zealand bacterial communities at the same order of magnitude of those found between the New Zealand streams, suggesting community composition was highly influenced by micro-environmental conditions within the mats. They also identified a strong and significant distance decay relationship (defined as an increasing dissimilarity in community composition with increasing geographic distance) between bacterial communities from New Zealand streams. However, the New Zealand samples were only taken at one site in each stream limiting knowledge on within-stream variability. Bouma-Gregson et al. (2018b) found an association with bacterial community composition within Microcoleus-dominated mats from Eel River (California, United States) and presence of Microcoleus strains capable of producing anatoxin-a (ATX). The authors also showed that many of the bacteria within the mats had metabolic capacities, such as oxygenic and anoxygenic photosynthesis, carbon respiration, sulfur compound oxidation and urea breakdown, which may benefit Microcoleus through the internal cycling of nutrients such as organic nitrogen.

      The overarching goal of this study was to investigate the bacterial community assemblages and anatoxin variability in M. autumnalis-dominated mats collected across different spatial scales. Five samples were collected from six sites of two streams in New Zealand and a suite of physical and chemical parameters were sampled or measured at each site. Bacterial communities were characterized using 16S rRNA metabarcoding and anatoxin content was determined using liquid chromatography-mass spectrometry (LC-MS). Droplet digital polymerase chain reaction (PCR) was used to determine anaC copies number which allowed anatoxin quotas (amount of anatoxin per cell) to be determined. We explored three questions; (1) how does the structure of bacterial communities in M. autumnalis-dominated mats vary within sites, within streams and between streams, (2) are anatoxin quotas as spatially variable as total anatoxin content in M. autumnalis-dominated mats across varying spatial scales, and (3) are there relationships between the structure of associated bacterial communities, anatoxin content and quota and stream water column physicochemical factors.

      Materials and Methods Sampling Sites

      Site surveys and M. autumnalis-dominated mat sample collection were undertaken at six sites (for geographic coordinates see Supplementary Material 1) along a 23 km longitudinal gradient of the Hutt River (Wellington, New Zealand) on November 28, 2017 and six sites over 9.5 km in the lower Cardrona River (Otago, New Zealand) on April 10–11, 2017 (Figure 1).

      Locations of the six study sites along the Hutt River (Wellington) and Cardrona River (Otago), New Zealand. Arrows indicate the direction of river flow.

      The Hutt River is a 56-km long aggrading stream (a stream that is actively elevating its bed by sediment deposition) with a catchment area of 655 km2. The upper reaches of the Hutt River catchment are principally native forest and scrublands. Land use in the lower reaches are comprised of low-intensity agricultural, commercial and residential areas.

      The Cardrona River has a total catchment area of 337 km2, it flows 40 km starting in the steep Cardrona Valley before, draining into the Clutha River. The Cardrona catchment primarily consists of agricultural grasslands ranging from tussock grasslands to high producing pasture grasslands in the upper and lower catchment, respectively. There is a heavy demand for water abstraction in the lower Cardrona River catchment, contributing to the lower stream drying in most years during summer. Upstream of site 2, the Cardrona River goes through a losing reach in which surface water is lost to groundwater with the stream only emerging ca. 400 m upstream of the sampling site. Approximately 1 km downstream the stream goes through another losing reach and was completely dry until emerging from the gravel approximately 200 m upstream of site 4. Site 1 and 6 were on small side stems of the main stream.

      Site Surveys and Sample Collection

      At each site, percentage coverage of M. autumnalis-dominated mats was visually assessed in five fields of view, equal distances apart, along three transects perpendicular to the bank edge using a bathyscope (Model 0800, Nuova Rade, Italy). At randomly chosen points along the transects, a small section (ca. 2 cm diameter) of a M. autumnalis-dominated mat was carefully removed using sterile tweezers, placed in a sterile tube (2 mL), and frozen (−20°C) for later anatoxin and DNA extractions. All mats were well developed and in their expansion phase (McAllister et al., 2016). This resulted in a total of five samples per site for all sites except for Hutt River site 4 (n = 4). A small subsample (0.5 g wet weight) of mat from each site was preserved in Lugol’s and water (5 mL). In the laboratory an inverted Olympus CKX-41 was used to confirm M. autumnalis-was the dominant cyanobacteria in the samples.

      Dissolved oxygen, conductivity (25°C reference temperature), pH and temperatures were measured at each site in the middle of the streams using a handheld water quality sonde (YSI Pro Plus, YSI Inc., Yellow Springs, OH, United States). The water velocity at each site was measured using a hand-held velocity meter (Marsh-McBirney, HACH, Loveland, CO, United States).

      Water samples (ca. 0.5–1 L) were collected adjacent to the most downstream mat at each site, syringe-filtered (Whatman GF/C, Whatman, Maidstone United Kingdom; ca. 1.6 μm pore size) and the filtrate frozen (−20°C) for dissolved nutrient analysis. The concentrations of nitrate (NO3-N), nitrite (NO2-N), ammoniacal-N (NH4-N) and dissolved reactive phosphorous (DRP) were determined at Hill Laboratories (Hamilton, New Zealand). Reporting limits for nutrient analyses were 0.001 mg L–1 (DRP, NO3-N, NO2-N, NH4-N) for Hutt River samples and 0.002 (NO2-N, NO3-N), 0.01 (NH4-N), and 0.004 mg L–1 (DRP) for the Cardrona River samples. Analysis was undertaken on a Lachat Quikchem® flow injection analyzer (FIA + 8000 Series, Zellweger Analytics, Inc.) using APHA (2012) 4500 methods. Dissolved inorganic nitrogen (DIN) was determined by summing NO3-N, NO2-N, and NH4-N.

      Anatoxin Extraction and Analysis

      Frozen M. autumnalis mat samples were defrosted, transferred to 20 mL glass vials and lyophilized (Gamma 1-16 LSC freeze-drier; Martin Christ, Germany). Lyophilized material was ground to a fine powder with a sterile metal spatula and two aliquots (for anatoxin and DNA extraction) were weighed (ca. 10–30 mg).

      The aliquots for anatoxin extraction were suspended in 1 mL of 0.1% formic acid (made up in Milli-Q water; Merck, Kenilworth, NJ, United States). The samples were vortexed (30 s), frozen (−20°C) and thawed in a sonicator bath (30 min). This was repeated two more times before the anatoxin extracts were centrifuged (12,000 × g, 5 min) and diluted 1/20 in 0.1% formic acid and frozen (−20°C). The supernatants were analyzed for dihydroanatoxin-a (dhATX), dihydrohomoanatoxin-a (dhHTX), ATX, and homoanatoxin-a (HTX) using LC-MS/MS as described in Wood et al. (2017b). Anatoxin concentrations were determined using an external standard curve constructed using dilutions of a certified reference material for ATX (National Research Council, Canada; 0.5–20 ng mL–1 in 0.1% formic acid). The ATX calibration curve was used to quantify HTX, dhATX, and dhHATX, using a relative response factor of 1. Data were converted to μg g–1 (equivalent to mg kg–1) by dividing the dilution-adjusted LC-MS results (in ng mL–1) by the weight of lyophilized starting material (in mg). Anatoxin quotas (in pg cell–1) were calculated by dividing the anatoxin concentration (in μg g–1; the sum of the four anatoxin congeners) by the concentration of anaC gene copies (in copies g–1; see below) and multiplying by 1,000 (to convert μg cell–1 to pg cell–1).

      DNA Extraction and <italic>anaC</italic> Gene Quantification

      The DNA aliquot was weighed directly into the first tube of the DNeasy PowerSoil® DNA Isolation Kit (QIAGEN, United States) and the extraction was performed following the manufacturer’s protocol. The quantity (>0.5 ng μL–1) of the extracted DNA was measured using a NanoPhotometer (Implen, Germany).

      Absolute concentrations of the anaC gene were measured using a BioRad QX200 Droplet Digital PCR system and M. autumnalis specific primers and a probe (Phor-AnaC-F5 5′-ACTAACCGAATCACTTCCACTT-3′, Phor-AnaC-R5 5′-CTC ACCCACCTCACCTTTAG-3′, Phor-AnaC-P5 5′-TTCAGTATT AGCGCAGGCTTTGCC-3′; Kelly et al., 2018). The hydrolysis probe was dual-labeled with a 5′ 6-carboxyfluorescein (6-FAM) fluorescent tag and a 3′ Black Hole Quencher. Each ddPCR reaction included 450 nM of each primer and probe, 1 × BioRad ddPCR Supermix for probes, 1 μL DNA, and sterile water for a total reaction volume of 22 μL. The BioRad QX200 droplet generator partitioned each reaction mixture into nanodroplets by combining 20 μL of the reaction mixture with 70 μL of BioRad droplet oil. After processing, a nanodroplet volume of 40 μL was transferred to a PCR plate for amplification using the following cycling protocol; hold at 95°C for 10 min, 40 cycles of 94°C for 30 s, 60°C 1 min, and a final enzyme deactivation step at 98°C for 10 min. The plate was then analyzed on the QX200 instrument (BioRad). For each ddPCR plate, at least one negative control (RNA/DNA-free water Life Technologies) and one positive control (genomic DNA extracted from a sample known to contain anatoxin) were included. When inhibition was observed, or samples were too concentrated, these were diluted with RNA/DNA-free water (UltraPureTM, Life Technologies, CA, United States) and reanalyzed. The results were then converted to copies g–1 using the following formula; number of copies per μL × 22 μL (the initial volume of the PCR reaction) × 100 μL (the volume used to elute the DNA during extraction)/weight of lyophilized starting material (in g).

      Bacterial Community Metabarcoding Polymerase Chain Reaction and High-Throughput Sequencing

      Amplification of PCR products for high-throughput sequencing was performed using bacterial specific primers to amplify a ca. 400 base pair (bp) sequence in the V3–V4 hypervariable region of the 16S rRNA gene; forward S-D-Bact-0341-b-S-17: 5′-CCT ACG GGN GGC WGC AG-3′ and reverse S-D-Bact-0785-a-A-21: 5′-GAC TAC HVG GGT ATC TAA TCC-3′ (Herlemann et al., 2011; Klindworth et al., 2013) modified to include IlluminaTM adapters. Each 50 μL PCR reaction contained 25 μL of AmpliTaq Gold® 360 master mix (Life Technologies, Camarillo, CA, United States), 5 μL of GC enhancer, 2 μL of each forward and reverse primers (10 μM), 15 μL of RNA/DNA free water (Thermo Fisher Scientific, United States) and 1 μL of template DNA. The PCR thermal cycling conditions involved denaturation for 10 min at 95°C, 27 cycles of denaturation at 95°C for 30 s, annealing at 50°C for 30 s, 72°C for 45 s and a final extension at 72°C for 7 min. PCR products were visualized by electrophoresis on 1.5% agarose gel with Red Safe DNA Loading Dye (iNtRON Biotechnology Inc., Kyungki-Do, South Korea) to confirm that only one PCR product was amplified at c. 400 bp. PCR products were cleaned (Agencourt® AMPure® XP Kit; Beckman Coulter, CA, United States), quantified (Qubit® 20 Fluorometer, Invitrogen), diluted to 10 ng μL–1 and submitted to New Zealand Genomics Limited (University of Auckland) for high-throughput Illumina sequencing.

      Bioinformatics

      Bioinformatics analysis was performed in QIIME2 (Caporaso et al., 2010). Primers were first removed from sequences, and sequencing data denoised and quality filtered with the default values of the DADA2 algorithm (Version 2018.2.6; Callahan et al., 2016). Amplicon sequence variants (ASVs) were taxonomically assigned with the q2-feature-classifier of QIIME2 and the SILVA reference database (Version 132; Yilmaz et al., 2014). Reference sequences were extracted and trimmed using the aforementioned primers. This was followed by training a multinomial Naïve Bayes classifier on the trimmed reference database and classifying representative sequences using the default values of QIIME2 feature-classifier. The raw FASTQ sequence data has been deposited in the Short Read Archive (NCBI project number PRJNA578643).

      Data Analysis

      Anatoxin quota, anaC copy number and total anatoxin concentrations were log-transformed prior to analysis in order to normalize their data distributions. One-way analyses of variance (ANOVAs) with Tukey Honest Significant Differences (HSD) pairwise comparisons were used to compare anatoxin concentrations, anatoxin quota, and anaC copy numbers between sites along either the Hutt or Cardrona rivers. Sites were assumed to be independent for all analyses. Welch Two Sample t-tests were used to compare toxin variables between the Hutt and Cardrona rivers. Statistical analyses were conducted using R software (R Studio Team, 2015; R Core Team, 2016) with ggplot2 (Wickham, 2016) and heplot (Fox et al., 2018).

      Prior to all analyses on bacterial community data, unassigned sequences, chloroplast, mitochondrial and cyanobacterial reads were removed from sequence data. ASVs present in method controls (n = 30) were removed from the dataset. Percentage stacked bar graphs were used to visualize the abundance of taxa in each sample using data prior to rarefaction. To assure adequate sequencing depth and normalize data information among samples, rarefaction curves were produced with the vegan package (Oksanen et al., 2018) using R software (Supplementary Material 2; R Studio Team, 2015; R Core Team, 2016). For alpha and beta diversity analysis sequencing depth was rarefied at 2,500 reads per sample through random subsampling with no replacement. Three samples with <2,500 reads were removed during rarefaction and excluded from further analysis (Hutt River 2B, 3E, and Cardrona River 4C).

      Beta-diversity assessment between streams was performed by transforming the rarefied data from both streams with a square root and Wisconsin standardization, followed by the production of a Bray-Curtis dissimilarity matrix. Based on this matrix, a non-metric multidimensional scaling (NMDS) was created to visualize the ordination of samples. Beta-diversity assessment (as previously described) was also carried out separately on Hutt River and Cardrona River data to determine finer-scale patterns. Statistical analyses were conducted using R software (R Studio Team, 2015; R Core Team, 2016) with ggplot2 (Wickham, 2016), phyloseq (McMurdie and Holmes, 2013), and vegan (Oksanen et al., 2018) packages.

      The multivariate differences in bacterial community composition between the two streams as well as between the sites in each stream was assessed using distance-based permutational analysis (PERMANOVA; Anderson et al., 2008) followed by pairwise tests. Analyses were based on Bray-Curtis similarities using square-root transformed data that had been rarefied as described above. This was conducted using PRIMER 7 with PERMANOVA + add-on (PRIMER-E Ltd., Plymouth, United Kingdom).

      The envfit function (vegan package) was used to assess the correlation between the ordination of the bacterial communities and log-transformed environmental variables as well as anatoxin concentration and quota. Significant variables (p < 0.05) were overlaid on the nMDS (described above) as vectors with the length of arrow proportional to the correlation strength and the arrow direction indicating the relationship between bacterial community and environmental variable.

      Results <italic>Microcoleus autumnalis</italic>-Dominated Mat Cover

      Microcoleus autumnalis-dominated mat cover was greater than 20% at all sites except for Hutt River site 3, where cover was 5% and Cardrona site 6, where cover was 10%. The highest M. autumnalis-dominated mat cover in both streams, 75 and 80%, was recorded for Hutt River sites 1 and 2, and Cardrona River site 4, respectively (Supplementary Material 3).

      Environmental Factors

      Conductivity varied between all sites, with the lowest conductivity in the Hutt River measured at site 2 (96 μS cm–1) and the highest at site 5 (116 μS cm–1). In the Cardrona River, conductivity tended to increase from the upstream to downstream sites, from 95 μS cm–1 at site 1 and 131 μS cm–1 at site 6 (Supplementary Material 3).

      DIN concentrations in the Hutt River were higher at the lower three sites, specifically site 4 where DIN measured 0.19 mg L–1. Upstream Hutt River sites tended to have less DIN, with the lowest concentrations measured at site 2 and 3 (0.097 and 0.071 mg L–1). There was more variability in DIN concentrations in the Cardrona River, with lower concentrations (sites 2 and 3; 0.056 mg L–1) and much higher concentrations, 0.36 and 0.42 mg L–1 recorded at sites 4 and 5, respectively (Supplementary Material 3).

      In the Hutt River, water column DRP tended to decrease from upstream (c. 0.025 mg L–1) to the most downstream site (site 6; <0.001 mg L–1). Except for Cardrona River site 6, where DRP was 0.015 mg L–1, DRP was generally below <0.01 mg L–1 (Supplementary Material 3).

      In the Hutt River, water temperature increased from upstream to downstream sites; (21.1 and 23.6°C at sites 1 and 6, respectively), except for site 3 where temperature was 20.8°C. Cardona River had substantially lower water temperature (ranging from 14.9 to 11.5°C at sites 3 and 5, respectively; Supplementary Material 3).

      Point water velocity ranged from 0.42 m s–1 (site 4) to 0.16 m s–1 (site 5) at Cardrona River sites. Higher velocities were recorded in the Hutt River, ranging from 0.15 m s–1 (site 6) to 0.83 m s–1 (site 3; Supplementary Material 3).

      Anatoxin Concentrations, Quotas and Congeners in <italic>Microcoleus autumnalis</italic>-Dominated Mats

      Median total anatoxins concentrations in the Hutt and Cardrona rivers were 7.56 and 21.95 mg kg–1 dw, respectively. Despite this difference in medians, there were no significant differences between total anatoxin concentrations in the Hutt and Cardrona rivers (t(54) = 1.17, p = 0.24) due to high within-stream and within site variability. In general, median total anatoxin concentration (the sum of the four structural congeners) per site increased from upstream to downstream, in both the Hutt and Cardrona rivers, with Cardrona site 1 being an exception to this trend (Figures 2A,B). There were statistically significant differences in total anatoxin concentrations between sites for both the Hutt River (F5,23 = 39.37, p < 0.001) and Cardrona River (F5,24 = 22.08, p < 0.001), with the highest median concentrations, 710 mg kg–1 dried weight (dw) and 778 mg kg–1 dw, in the Hutt and Cardrona rivers, respectively, measured at Hutt site 6 and Cardrona site 6 (Figures 2A,B). Relatively high anatoxin concentrations were also measured at Cardrona site 1 (median: 647 mg kg–1 dw). The highest maximum total anatoxin concentrations in each stream were measured in samples from Hutt River site 4 (2,116 mg kg–1 dw) and Cardrona River site 6 (1,017 kg–1 dw; Figures 2A,B).

      (A,B) Total anatoxin (i.e., the sum of the four anatoxin congeners); (C,D) concentrations of anaC copy numbers; and (E,F) anatoxin quotas (amount of toxin per toxic cell) from Microcoleus autumnalis-dominated mats collected from sites along the Hutt River (Wellington, New Zealand) or the Cardrona River (Otago, New Zealand). Sites are arranged along x-axis in an upstream (left) to downstream (right) gradient. n = 5 for all sites except for Hutt River site 4 (n = 4). Note that y-axes are different log scales. Solid black lines indicate the median, box shows 1st and 3rd quartiles, whiskers extend to the last data point within 1.5 times the inter-quartile range. Open circles are outliers beyond this range. Shared letters indicate no significant difference (p > 0.05) between sites (One-way ANOVA, Tukey HSD). Dw, dry weight; pg, picograms.

      Comparison of the anaC gene copy numbers between Hutt and Cardrona rivers revealed a significant difference (t(55) = -2.24, p = 0.03), with higher copy numbers in Hutt River samples (median = 135) in comparison to Cardrona River (median = 68.0). In the Hutt River, anaC gene copy numbers differed significantly between sites (F5,23 = 39.72, p < 0.001) with downstream sites 4, 5, and 6 having significantly higher (at least 20- and 200-fold higher, respectively) copy numbers than upstream sites (Tukey HSD, p < 0.001; Figure 2C). Similarly, in the Cardrona River there was significant variation in anaC gene copies between sites (F5,24 = 8.80, p < 0.001); however, this trend was not longitudinal with sites 1 and 6 having at least 15- and 10-fold higher copy numbers than other sites, respectively (Figure 2D).

      Between the Hutt and Cardrona rivers, there were significant differences in anatoxin quotas (t(51) = 4.50, p < 0.001), with quotas higher in the Cardrona River (median = 0.40 pg cell–1) compared to the Hutt River (median = 0.054 pg cell–1). Anatoxins quotas differed significantly among Hutt River sites (F5,23 = 20.76, p < 0.001) and Cardrona River sites (F5,24 = 8.12, p < 0.001). In the Hutt River, quotas were significantly higher at sites 3 (median 0.44 pg cell–1) and site 4 (median 0.58 pg cell–1) in comparison to other sites, specifically site 2 (median 0.002 pg cell–1; Tukey HSD, p < 0.01; Figure 2E). In the Cardrona River, site 5 and 6 had the highest median quotas (1.82 and 1.16 pg cell–1, respectively) and the lowest median quota was measured in site 3 samples (0.04 pg cell–1). The highest quota values overall were measured in Cardrona River samples (6.45 and 11.72 pg cell–1, sites 4 and 5, respectively; Figure 2F).

      Dihydro-ATX was the most abundant congener detected in the M. autumnalis-dominated samples at Cardrona River sites (58–87%) and at all Hutt River sites (>47%) except for site 1, where HTX was most abundant (58%; Figures 3A,B). In the Hutt River, HTX was the second most abundant anatoxin congener at all sites, contributing 6–58% (Figure 3A). In contrast, dhHTX was the second most abundant congener in the Cardrona River, contributing between 7–20% of total toxins at each site (Figure 3B). ATX contributed less than 1% abundance to the congener profiles of samples from the Cardrona and Hutt rivers (Figures 3A,B).

      Stacked bar graph showing the average percentage of each anatoxin congener in Microcoleus autumnalis-dominated mat samples collected from (A) Hutt River (Wellington, New Zealand) and (B) Cardrona River (Otago, New Zealand). ATX, anatoxin-a; HTX, homoanatoxin-a; dhATX, dihydroanatoxin-a; dhHTX, dihydrohomoanatoxin-a. n = 5 for all sites except for Hutt River site 4 (n = 4). Sites are arranged along x-axis in an upstream (left) to downstream (right) gradient.

      Bacterial Communities in <italic>Microcoleus autumnalis</italic>-Dominated Mats

      Analysis of the 16S rRNA gene sequences showed that all the samples were dominated by cyanobacterial reads (mean ± SD: 76.0 ± 18.3% of the total reads per site) with Phormidiaceae contributing 98.7 ± 2.6% of total cyanobacterial reads. The remaining 1.3% cyanobacterial reads were primarily Pseudanabaenales followed by Leptolyngbyales and to a lesser degree, Synechococcales and Caenarcaniphilales (Melainabacteria).

      Following the removal of cyanobacterial sequences and post-rarefaction, a total of 1372 and 1098 ASVs were present across all Hutt and Cardrona river samples, respectively (Figure 4A). In rarefied data, only 11.2% of all ASVs were shared among the Hutt and Cardrona samples with the remaining 88.8% unique to either stream (Figure 4A). Upon examination of the shared ASVs (Figure 4B), patterns in bacterial composition were consistent with those seen when the total bacterial communities were analyzed (Figure 5).

      (A) Venn diagram of bacterial amplicon sequence variants (ASVs) from Microcoleus autumnalis-dominated mats from the Hutt and Cardrona rivers, New Zealand. ASVs are based on 16S rRNA gene fragment sequences. Number of ASVs in each river and shared ASVs between rivers are shown excluding chloroplast and cyanobacterial reads for post rarefaction of sample sizes to 2500 reads per sample. (B) Stacked bar graph showing the relative abundance of the shared ASVs from the bacterial 16S rRNA sequences in samples from five Hutt River (Wellington, New Zealand) sites, and six Cardrona River (Otago, New Zealand) sites. The most abundant phyla (contributing >1% of total abundance in each river) are stated. “Other” composed of rare phyla (<1%).

      Stacked bar graph showing the relative abundance of bacterial 16S rRNA sequences at (A) phyla, (B) class, (C) order taxonomy levels in Microcoleus autumnalis-dominated samples from five Hutt River (Wellington, New Zealand) sites, and six Cardrona River (Otago, New Zealand) sites. Cyanobacterial reads were excluded, and the most abundant phyla (contributing >1% of total abundance in each river) are stated. “Other” composed of rare phyla (<1%). Sites are arranged along x-axis in an upstream (left) to downstream (right) gradient.

      When examining the total bacterial communities in both streams but following the exclusion of all cyanobacterial reads, the most abundant bacterial phyla present was Bacteroidetes, contributing 69.9 ± 11.2% and 64.7 ± 17.8% of the total abundance of the bacterial community in samples from the Hutt and Cardrona rivers, respectively (Figure 5A). Among Bacteroidetes, Bacteroidia dominated the class level, composing 69.4 ± 11.4% (of total abundance) in the Hutt River and 64.4 ± 17.8% in the Cardrona River (Figure 5B). The top three orders of Bacteroidetes in order of abundance in the Hutt River were Sphingobacteriales (26.5 ± 10.1% of total abundance), Cytophagales (20.7 ± 9.9%), and Chitinophagales (15.9 ± 4.9%) with Flavobacteriales present at a lower abundance (4.7 ± 3.2%; (Figure 5C). Pedobacter (Sphingobacteriales; Bacteroidetes) was the most abundant genus (20.9 ± 11.6%) in the Hutt River. Similar bacterial orders were dominant in the Cardrona River; however, with some differences in respective abundances with successive dominant orders from Cytophagales (20.9 ± 12.2%), Chitinophagales (20.9 ± 7.6%), and Sphingobacteriales (14.1 ± 5.9%; Figure 5C). Flavobacteriales were present at lower abundance (4.5 ± 2.7%) in the Cardrona River. The most abundant genera in the Cardrona River differed to the Hutt River with Ferruginibacter (8.2 ± 4.5%) and Sediminibacterium (7.8 ± 5.1%) from the Chitinophagales order dominating.

      Proteobacteria was the second most abundant phylum, making up 22.5 ± 8.5% of the Hutt River and 29.4 ± 14.2% of the Cardrona River total bacterial community composition. Alphaproteobacteria was dominant within Proteobacteria, contributing 20.1 ± 8.3% of total abundance in the Hutt River whereas although present, Deltaproteobacteria contributed less (2.1 ± 2%) to the total abundance (Figure 5B). Within Alphaproteobacteria, most of the sequences in the Hutt River samples were classified in the orders of Acetobacterales (5.3 ± 3.2% of total abundance), Sphingomonadales (5.0 ± 2.7%), and Rhodobacterales (3.3 ± 2.3%; Figure 5C). A similar pattern was observed in the Cardrona River with Alphaproteobacteria the dominant Proteobacteria, composing 27.9 ± 14.1% of total abundance (Figure 5C). In the Cardrona, Rhodobacterales and Sphingomonadales were the most abundant Alphaproteobacteria, contributing 10.2 ± 6.6% and 7.6 ± 5.2% of total abundance, respectively.

      Multidimensional scaling based on Bray-Curtis similarities between bacterial communities (using ASVs), showed a clear separation between the two streams (Figure 6A). Analysis of overall bacterial community structure using PERMANOVA revealed a significant difference between the Hutt and Cardrona rivers (pseudo-F = 18.4, p < 0.001). Due to this high degree of separation, the relationship between bacterial communities and the environmental variables and anatoxin concentration was also compared for each stream separately (Figures 6B,C).

      Non-metric, multi-dimensional scaling ordination (confidence level = 0.95) based on Bray-Curtis similarities with 999 permutations of bacterial communities in Microcoleus autumnalis-dominated mat samples from; (A) the Hutt and Cardrona rivers, (B) the Cardrona River, and (C) the Hutt River at the individual sequence level. Relationships with dissolved reactive phosphorus (DRP), conductivity, velocity, dissolved inorganic nitrogen (DIN), temperature, and total anatoxin concentrations (log transformed) are indicated by vectors. Only significant relationships (p < 0.05) are shown.

      Bacterial communities from the Hutt River sites 1, 2, and 3 (upstream sites) formed a distinct cluster, separate to the mid-stream site 4 and the downstream sites 5 and 6, which also grouped together (Figure 6C). PERMANOVA analysis showed that the bacterial communities were significantly different between Hutt River sites (pseudo-F = 5.14, p < 0.001). Pairwise analysis revealed that all sites were significantly different, with the only exception being sites 2 and 3 (Supplementary Material 4). In the Cardrona River, except for an overlap between two sites (sites 2 and 5), each site clustered distinctly with significant differences between sites (PERMANOVA; pseudo-F = 4.78, p < 0.001; Figure 6B). Furthermore, pairwise tests showed that despite this overlap, the bacterial communities at all sites were significantly different from each other (Supplementary Material 5).

      Relationship Between Water Chemistry, Total Anatoxins and Quota’s and Bacterial Communities

      When comparing between streams, the bacterial community composition in the Hutt River was correlated with the higher velocity and temperature (Figure 6A). In comparison, higher DRP, anatoxin quota and concentration were associated with bacterial community composition in the Cardrona River (Figure 6A). In both streams, conductivity was also significantly associated with the distribution of bacterial communities between sites.

      The environmental variables that were strongly associated with bacterial community compositions in the Cardrona River were conductivity and DRP, with higher levels correlated with site 4 and 6, respectively (Figure 6B). There was also a relationship between anatoxin concentration/quotas and bacterial communities in the Cardrona River; however, this relationship was not specific to one site.

      The composition of clustered bacterial communities identified at the upper sites 1, 2, and 3 in the Hutt River were associated with the higher DRP and water velocity found at these sites in comparison to downstream sites (Figure 6C). Higher DIN, conductivity and temperature were correlated with the bacterial community composition at Hutt River site 6 (most downstream site). Additionally, higher total anatoxin concentrations as well as anatoxin quota were strongly associated with the composition of the bacterial communities at Hutt River sites 4 and 6.

      Discussion

      Using 16S rRNA metabarcoding we characterized diverse bacterial assemblages in M. autumnalis-dominated mats from two large fluvial systems (the Hutt and Cardrona rivers) in New Zealand. This data, in concert with previous studies (Brasell et al., 2015; Bouma-Gregson et al., 2018b; Echenique-Subiabre et al., 2018) clearly demonstrate that Microcoleus-dominated mats are complex microbial communities. Although there were fine-scale differences in bacterial community composition between streams, sites and samples when analyzed at the ASV taxonomic level; at higher taxonomic ranks (i.e., phyla) both the Hutt and Cardrona river bacterial communities were largely composed of Proteobacteria and Bacteroidetes. This is similar to that described by Brasell et al. (2015) who worked on a single site in the Hutt River, and Echenique-Subiabre et al. (2018) who investigated Microcoleus-dominated mats in the Tarn River in France and eight New Zealand rivers and streams. In contrast to these studies, the present research showed that Bacteroidetes dominated whereas, both Brasell et al. (2015) and Echenique-Subiabre et al. (2018) identified Proteobacteria as the dominant phylum. Furthermore, both Betaproteobacteria and Alphaproteobacteria were the main and largely equal contributors to Proteobacterial abundance in the studies of Brasell et al. (2015) and Bouma-Gregson et al. (2018b), differing with the almost exclusive presence of Alphaproteobacteria in the present study. As a result of this, Burkholderiales (Betaproteobacteria), although present in samples, did not contribute the bacterial composition to the same degree that was identified by both Brasell et al. (2015) and Echenique-Subiabre et al. (2018). Alphaproteobacteria orders Sphingomonadales and Rhodobacterales were prominent across all studies (Brasell et al., 2015; Echenique-Subiabre et al., 2018). Bouma-Gregson et al. (2018b) worked on Microcoleus-dominated mats in the Eel River, California and the abundant phyla described in their manuscript are more aligned to the present study with Bacteroidetes dominant. Differences between studies could also be partly due to biases that are introduced by employing different high-throughput sequencing and bioinformatic methodologies; for example, the selection of DNA metabarcoding region and primers, sequencing platforms, the preparation of sequence libraries including the reference database used for taxonomic assignment and subsequent clustering methods (eg., Balvočiūtė and Huson, 2017; Abdelfattah et al., 2018; Porter and Hajibabaei, 2018; Bukin et al., 2019; Zinger et al., 2019). A limitation of the current study is the low rarefaction threshold, which may reduce the overall coverage of biodiversity compared to what has been identified in other studies.

      There were distinct bacterial communities identified at a range of spatial scales with significant differences identified between streams and sites. The between stream differences were apparent despite the same organism, Microcoleus, dominating the mat community in both streams. In addition to environmental factors (discussed below) there are several considerations regarding the influence of dominant taxa in the mat assemblage. The two streams were sampled at different times of the year; early summer and late autumn (Hutt and Cardrona rivers, respectively). It is possible the residence period of the mats prior to sample collection differed between each river due to differences in seasonal growth and flushing flows. Brasell et al. (2015) described distinct phases of succession in microbial communities associated with Microcoleus, highlighting that mat age is an important consideration that may explain differences between stream assemblages. There is also a potential that Microcoleus is genetically different between streams (this study only examined the 16S rRNA gene). Genetic variation in functional genes can result in biological consequences, for example variability in anatoxin concentrations among Microcoleus-dominated mats (Wood and Puddick, 2017), that may in turn influence the associated microbial community. Additionally, there may be other co-inhabiting organisms in the mats (i.e., diatoms or other eukaryotic algae) that differ between streams and sites and could influence co-inhabiting bacteria. Similar patterns found by Echenique-Subiabre et al. (2018) showed greater dissimilarity between New Zealand streams, compared to dissimilarity between New Zealand streams and a French stream, which led them to suggest that factors within mat habitats were more likely to regulate associated microbial community in comparison to external factors.

      At a finer spatial scale, some sites exhibited considerable within-site variability that was sometimes greater than between-site variability. Additionally, when there was considerable overlap in communities between sites (i.e., Cardrona River sites 2 and 5), these patterns did not consistently follow spatial patterns such as a predicted higher similarity between sites situated next to each other. Stream water potentially provides an inoculum for microbial communities; however, the results of this study indicate this may not be the key inoculum. This was especially apparent in the Cardrona River between sites 2 and 3, which are very close to each other and connected by flowing water but isolated from all other sites (due to the stream being dry, i.e., losing reach). Despite this, the communities in the M. autumnalis-dominated mats at these sites were significantly different. Cardrona sites 4, 5, and 6 were all connected by flowing water, but there were significant differences in bacterial community composition between these sites. This suggests other environmental parameters either external or within mats are playing a significant role in driving community structure.

      It is well known that environmental conditions can influence the bacterial community composition of stream bacteria (e.g., Fierer et al., 2007; Zeglin, 2015). The data in this study indicates that temperature, water chemistry, concentration of DRP as well as conductivity, and water velocity may play a role in structuring bacterial communities. However, caution is advised when relating water column measurements to within mat communities. Conditions in M. autumnalis-dominated mats can be markedly different to those of the overlying stream water. For example, Wood et al. (2015) showed that water within M. autumnalis-dominated mats can have up to 320-fold higher DRP concentrations compared to surrounding stream water as a result of desorption of phosphorus bound to sediment within the mat matrix as pH and dissolved oxygen changed markedly over diurnal cycles and were very different to water column values. It is likely that micro-environmental conditions, such as these, have a substantial effect on bacterial community structure, but that these organisms also play an important role in creating and maintaining these conditions. Bouma-Gregson et al. (2018b) suggested, based on metagenomic analyses, that the microbial community contributes to many processes such as phosphorus mineralization in mats through the scavenging of phosphorus using acid phosphatases during periods of pH elevation. Nitrogen cycling within mats, with the degradation of urea as a nitrogenous waste product recycled by microbial organisms, is also suggested to decouple nitrogen requirements from overlying contributions from the external system (Bouma-Gregson et al., 2018b); a state that is posited to occur when mats are mature and established. Further studies, which use techniques that enable conditions within M. autumnalis-dominated mats to be characterized are required to fully understand the relationship between within-mat conditions and bacterial community composition. In addition to exploring structural differences among communities, enhanced knowledge on their functions (i.e., metatranscriptomics) would provide valuable insights as to whether communities are involved in similar processes in the mat environments, despite differences in community composition.

      Statistical analysis identified a relationship between bacterial community structure in M. autumnalis-dominated mats and total anatoxin concentration as well as anatoxin quota. Bouma-Gregson et al. (2018b) also identified a difference between the bacterial communities from Microcoleus-dominated mat samples with and without the anatoxin-a gene operon. Although it was possible that bacteria present were influencing selection for toxin-producing cyanobacterial strains, Bouma-Gregson et al. (2018b) considered it more likely that cyanobacterial-produced compounds, such as anatoxin-a, were partly responsible for shaping the microbial community. Harland et al. (2013) and Harland et al. (2014) hypothesized that cyanotoxins may be involved in intra- or inter-species interactions, specifically chemotaxic signaling either between M. autumnalis cells and/or other micro-organisms. Their laboratory-based studies showed that up-regulation in anatoxin production occurred early in the growth phase, which could indicate a role for the toxin in the colonization stage of cyanobacterial mats.

      The anatoxin concentration data collected in the present study conformed to patterns observed in other studies, with considerable spatial variability between mats at each site and between sites (Wood et al., 2010, 2017a; McAllister et al., 2018). In general, the highest median anatoxin concentration values (ca. 700 mg kg–1) for each stream were markedly higher than most of those reported for other New Zealand streams (McAllister et al., 2016). The relative abundance of anatoxin congeners observed in the samples collected during the present study generally contained dhATX as the dominant congener followed by HTX or dhHTX. This was consistent with results from other New Zealand streams (McAllister et al., 2016). This highlights the importance of analyzing for all anatoxin congeners (not solely ATX) when undertaking research on M. autumnalis and when assessing the human and animal health risks.

      Although numerous studies have shown high variability in anatoxin content of Microcoleus-dominated mats, few studies have explored differences in anatoxin quotas (Wood and Puddick, 2017; Kelly et al., 2018). In general, converting the data to anatoxin quotas reduced the variability among the mat samples in comparison to total toxin concentrations. This corroborates the findings of Wood and Puddick (2017) who demonstrated that the abundance of toxic genotypes within M. autumnalis-dominated mats is likely to be the main cause of variability in total anatoxin concentrations among mats. Although variability was reduced by the conversion of total toxin to anatoxin quota, some differences remained. These could be due to an up- or down-regulation in anatoxin production caused by biotic or abiotic factors (discussed below), however, given that no consistent patterns were observed this is more likely due to differences in the amount of toxin produced by different strains of toxic M. autumnalis cohabitating in the mats. Wood et al. (2012) isolated and cultured 30 strains of M. autumnalis from four small sections of benthic mats. Among the toxic strains there were 100-fold differences in the amount of anatoxin produced.

      Although environmental factors regulating anatoxin production were not directly investigated in this study, previous research by Heath et al. (2011) suggested that toxin-producing strains in the Hutt River “out-competed” non-toxic M. autumnalis-dominated strains at temperatures above 15°C, but subsequent studies in other streams have not observed this pattern (Wood et al., 2017a). Heath et al. (2014) investigated the effects of nitrogen and phosphorus on the growth of two M. autumnalis strains (toxic and non-toxic) in a laboratory-based study. Cell concentrations and maximum growth rates were higher for the non-anatoxin-producing strain regardless of treatments, suggesting that specific genotypes may have environmental preferences. Experiments using a larger number of strains are required before definitive conclusions can be drawn about environmental preferences of toxic versus non-toxic M. autumnalis genotypes.

      Conclusion

      This study continues to highlight the complexity of the microbial communities within M. autumnalis-dominated mats. Some of these species will be involved in internal nutrient cycling processes, potentially explaining why establishing relationships between the abundance of M. autumnalis-dominated mats and environmental variables measured in stream water is complex or impossible. Further studies investigating the functions associated with these communities are recommended. The results of the present study corroborated previous research and show high variability in total anatoxin concentrations at varying spatial scales. When the data were converted to anatoxin quotas this variability was reduced, suggesting that the relative abundance of toxic genotypes within Microcoleus-dominated mats is a key driver of differences in total anatoxin concentrations among mats. Factors that regulate the relative abundance of toxic and non-toxic genotypes or anatoxin production remain unknown.

      Data Availability Statement

      The datasets generated for this study can be found in NCBI SRA accession PRJNA578643.

      Author Contributions

      SW, JP, and MH contributed to the conception and design of the study. SF, SW, JP, and KS acquired the samples and data. OL performed the statistical analysis. GT-L, OL, SF, and SW contributed to the analysis and interpretation of the data. GT-L and SW wrote the first draft of the manuscript. All authors contributed to the manuscript revision, read, and approved the submitted version.

      Conflict of Interest

      The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

      Funding. SW received funding from the National Institute of Water and Atmospheric Research Ltd. under the causes and effects of water quality degradation: eutrophication risk assessment program. The New Zealand Ministry of Environment provided funding and support to undertake this work.

      We thank Joseph Fraser, Adam Uytendaal, Dean Olsen (Otago Regional Council), Katie Brasell, Emily Martin (Greater Wellington Regional Council), Charlotte Tomlinson (Cawthron Foundation 2017/2018 Sir Theodore Rigg Scholar), and Laura Kelly (Victoria University) for sample collection and field assistance. Jacob Thomson-Laing, Charlotte Tomlinson, and Carrie Paige (Cawthron) assisted with the preparation of samples from this study. We thank Evan Harrison (Ministry for the Environment/Greater Wellington Regional Council) for his support and advice throughout the project.

      Supplementary Material

      The Supplementary Material for this article can be found online at: /articles/10.3389/fmicb.2020.00129/full#supplementary-material

      References Abdelfattah A. Malacrinò A. Wisniewski M. Cacciola S. O. Schena L. (2018). Metabarcoding: a powerful tool to investigate microbial communities and shape future plant protection strategies. Biol. Control 120 110. 10.1016/j.biocontrol.2017.07.009 Aboal M. Puig M. Á (2005). Intracellular and dissolved microcystin in reservoirs of the river Segura basin, Murcia, SE Spain. Toxicon 45 509518. 10.1016/j.toxicon.2004.12.012 15733573 Anderson M. Gorley R. N. Clarke K. (2008). PERMANOVA+ for Primer: Guide to Software and Statistical Methods. Plymouth: PRIMER-E. Balvočiūtė M. Huson D. H. (2017). SILVA, RDP, Greengenes, NCBI and OTT — how do these taxonomies compare? BMC Genomics 18:114. 10.1186/s12864-017-3501-3504 28361695 Berry M. A. Davis T. W. Cory R. M. Duhaime M. B. Johengen T. H. Kling G. W. (2017). Cyanobacterial harmful algal blooms are a biological disturbance to Western Lake Erie bacterial communities. Environ. Microbiol. 19 11491162. 10.1111/1462-2920.13640 28026093 Borges H. L. F. Branco L. H. Z. Martins M. D. Lima C. S. Barbosa P. T. Lira G. A. S. T. (2015). Cyanotoxin production and phylogeny of benthic cyanobacterial strains isolated from the northeast of Brazil. Harmf. Algae 43 4657. 10.1016/j.hal.2015.01.003 Bouma-Gregson K. Kudela R. M. Power M. E. (2018a). Widespread anatoxin-a detection in benthic cyanobacterial mats throughout a river network. PLoS One 13:e0197669. 10.1371/journal.pone.0197669 29775481 Bouma-Gregson K. Olm M. R. Probst A. J. Anantharaman K. Power M. E. Banfield J. F. (2018b). Microbial diversity and metabolic potential in cyanotoxin producing cyanobacterial mats throughout a river network. bioRxiv [Preprint] 10.1371/journal.pone.0197669 29775481 Brasell K. A. Heath M. W. Ryan K. G. Wood S. A. (2015). Successional change in microbial communities of benthic Phormidium-dominated biofilms. Microb. Ecol. 69 254266. 10.1007/s00248-014-0538-537 25467742 Bukin Y. S. Galachyants Y. P. Morozov I. V. Bukin S. V. Zakharenko A. S. Zemskaya T. I. (2019). The effect of 16S rRNA region choice on bacterial community metabarcoding results. Sci. Data 6:190007. 10.1038/sdata.2019.7 30720800 Cai H. Jiang H. Krumholz L. R. Yang Z. (2014). Bacterial community composition of size-fractioned aggregates within the phycosphere of cyanobacterial blooms in a eutrophic freshwater lake. PLoS One 9:e102879. 10.1371/journal.pone.0102879 25144467 Callahan B. J. McMurdie P. J. Rosen M. J. Han A. W. Johnson A. J. Holmes S. P. (2016). DADA2: high-resolution sample inference from Illumina amplicon data. Nat. Methods 13 581583. 10.1038/nmeth.3869 27214047 Caporaso J. G. Kuczynski J. Stombaugh J. Bittinger K. Bushman F. D. Costello E. K. (2010). QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7:335. 10.1038/nmeth.f.303 20383131 Chorus I. Bartram J. (1999). Toxic Cyanobacteria in Water: a Guide to Their Public Health Consequences, Monitoring and Management. Geneva: World Health Organization. Cirés S. Casero M. Quesada A. (2017). Toxicity at the edge of life: a review on cyanobacterial toxins from extreme environments. Mar. Drugs 15:233. 10.3390/md15070233 28737704 Codd G. A. (2000). Cyanobacterial toxins, the perception of water quality, and the prioritisation of eutrophication control. Ecol. Eng. 16 5160. 10.1016/S0925-8574(00)00089-86 Codd G. A. Bell S. G. Kaya K. Ward C. J. Beattie K. A. Metcalf J. S. (1999). Cyanobacterial toxins, exposure routes and human health. Eur. J. Phycol. 34 405415. 10.1080/09670269910001736462 Echenique-Subiabre I. Zancarini A. Heath M. W. Wood S. A. Quiblier C. Humbert J. F. (2018). Multiple processes acting from local to large geographical scales shape bacterial communities associated with Phormidium (cyanobacteria) biofilms in French and New Zealand rivers. Sci. Rep. 8:14416. 10.1038/s41598-018-32772-w 30258224 Faassen E. J. Harkema L. Begeman L. Lurling M. (2012). First report of (homo)anatoxin-a and dog neurotoxicosis after ingestion of benthic cyanobacteria in The Netherlands. Toxicon 60 378384. 10.1016/j.toxicon.2012.04.335 22534073 Fetscher A. E. Howard M. D. A. Stancheva R. Kudela R. M. Stein E. D. Sutula M. A. (2015). Wadeable streams as widespread sources of benthic cyanotoxins in California, USA. Harmf. Algae 49 105116. 10.1016/j.hal.2015.09.002 Fierer N. Morse J. L. Berthrong S. T. Bernhardt E. S. Jackson R. B. (2007). Environmental controls on the landscape-scale biogeography of stream bacterial communities. Ecology 88 21622173. 10.1890/06-1746.1 17918395 Foss A. J. Butt J. Aubel M. T. (2018). Benthic periphyton from Pennsylvania, USA is a source for both hepatotoxins (microcystins/nodularin) and neurotoxins (anatoxin-a/homoanatoxin-a). Toxicon 150 1316. 10.1016/j.toxicon.2018.05.002 29746979 Foss A. J. Butt J. Fuller S. Cieslik K. Aubel M. T. Wertz T. (2017). Nodularin from benthic freshwater periphyton and implications for trophic transfer. Toxicon 140 4559. 10.1016/j.toxicon.2017.10.023 29107081 Fox J. Friendly M. Monette G. (2018). Heplots: Visualizing Tests in Multivariate Linear Models. R package version 1.3–5. Gugger M. Lenoir S. Berger C. Ledreux A. Druart J.-C. Humbert J.-F. (2005). First report in a river in France of the benthic cyanobacterium Phormidium favosum producing anatoxin-a associated with dog neurotoxicosis. Toxicon 45 919928. 10.1016/j.toxicon.2005.02.031 15904687 Harland F. M. J. Wood S. A. Broady P. A. Gaw S. Williamson W. M. (2014). Polyphasic studies of cyanobacterial strains isolated from benthic freshwater mats in Canterbury, New Zealand. New Zealand J. Bot. 52 116135. 10.1080/0028825x.2013.846266 Harland F. M. J. Wood S. A. Moltchanova E. Williamson W. M. Gaw S. (2013). Phormidium autumnale growth and anatoxin- a production under iron and copper stress. Toxins 5 25042521. 10.3390/toxins5122504 24351714 Hart D. D. Biggs B. J. F. Nikora V. I. Flinders C. A. (2013). Flow effects on periphyton patches and their ecological consequences in a New Zealand river. Freshw. Biol. 58 15881602. 10.1111/fwb.12147 Heath M. W. Wood S. A. Barbieri R. F. Young R. G. Ryan K. G. (2014). Effects of nitrogen and phosphorus on anatoxin-a, homoanatoxin-a, dihydroanatoxin-a and dihydrohomoanatoxin-a production by Phormidium autumnale. Toxicon 92 179185. 10.1016/j.toxicon.2014.10.014 25449104 Heath M. W. Wood S. A. Ryan K. G. (2011). Spatial and temporal variability in Phormidium mats and associated anatoxin-a and homoanatoxin-a in two New Zealand rivers. Aquatic Microb. Ecol. 64 6979. 10.3354/ame01516 Herlemann D. P. R. Labrenz M. Jürgens K. Bertilsson S. Waniek J. J. Andersson A. F. (2011). Transitions in bacterial communities along the 2000 km salinity gradient of the Baltic Sea. ISME J. 5:1571. 10.1038/ismej.2011.41 24928874 Huisman J. Codd G. A. Paerl H. W. Ibelings B. W. Verspagen J. M. H. Visser P. M. (2018). Cyanobacterial blooms. Nat. Rev. Microbiol. 16 471483. 10.1038/s41579-018-0040-41 29946124 Kelly L. T. Wood S. A. McAllister T. G. Ryan K. G. (2018). Development and application of a quantitative PCR assay to assess genotype dynamics and anatoxin content in Microcoleus autumnalis-dominated mats. Toxins 10:431. 10.3390/toxins10110431 30373141 Klindworth A. Pruesse E. Schweer T. Peplies J. Quast C. Horn M. (2013). Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 41:e1. 10.1093/nar/gks808 22933715 Lee J. Lee S. Jiang X. (2017). Cyanobacterial toxins in freshwater and food: important sources of exposure to humans. Ann. Rev. Food Sci. Technol. 8 281304. 10.1146/annurev-food-030216-30116 28245155 Louati I. Pascault N. Debroas D. Bernard C. Humbert J.-F. Leloup J. (2015). Structural diversity of bacterial communities associated with bloom-forming freshwater cyanobacteria differs according to the cyanobacterial genus. PLoS One 10:e0140614. 10.1371/journal.pone.0140614 26579722 McAllister T. G. Wood S. A. Atalah J. Hawes I. (2018). Spatiotemporal dynamics of Phormidium cover and anatoxin concentrations in eight New Zealand rivers with contrasting nutrient and flow regimes. Sci. Total Environ. 612 7180. 10.1016/j.scitotenv.2017.08.085 28846906 McAllister T. G. Wood S. A. Hawes I. (2016). The rise of toxic benthic Phormidium proliferations: a review of their taxonomy, distribution, toxin content and factors regulating prevalence and increased severity. Harmf. Algae 55 282294. 10.1016/j.hal.2016.04.002 28073542 McMurdie P. J. Holmes S. (2013). Phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8:e61217. 10.1371/journal.pone.0061217 23630581 Mez K. Beattie K. A. Codd G. A. Hanselmann K. Hauser B. Naegeli H. (1997). Identification of a microcystin in benthic cyanobacteria linked to cattle deaths on alpine pastures in Switzerland. Eur. J. Phycol. 32 111117. 10.1080/09670269710001737029 Mohamed Z. A. El-Sharouny H. M. Ali W. S. M. (2006). Microcystin production in benthic mats of cyanobacteria in the Nile River and irrigation canals Egypt. Toxicon 47 584590. 10.1016/j.toxicon.2006.01.029 16564062 Mowe M. A. D. Mitrovic S. M. Lim R. P. Furey A. Yeo D. C. J. (2015). Tropical cyanobacterial blooms: a review of prevalence, problem taxa, toxins and influencing environmental factors. J. Limnol. 74 205224. 10.4081/jlimnol.2014.1005 Oksanen J. Blanchet F. G. Friendly M. Kindt R. Legendre P. McGlinn D. (2018). Vegan: Community Ecology Package. R Package Version 2.5-1. O’Neil J. M. Davis T. W. Burford M. A. Gobler C. J. (2012). The rise of harmful cyanobacteria blooms: the potential roles of eutrophication and climate change. Harmf. Algae 14 313334. 10.1016/j.hal.2011.10.027 Paerl H. W. Fulton R. S. III Moisander P. H. Dyble J. (2001). Harmful freshwater algal blooms, with an emphasis on cyanobacteria. Sci. World J. 1 76113. 10.1100/tsw.2001.16 12805693 Paerl H. W. Paul V. J. (2012). Climate change: links to global expansion of harmful cyanobacteria. Water Res. 46 13491363. 10.1016/j.watres.2011.08.002 21893330 Paerl H. W. Pinckney J. L. Steppe T. F. (2000). Cyanobacterial–bacterial mat consortia: examining the functional unit of microbial survival and growth in extreme environments. Environ. Microbiol. 2 1126. 10.1046/j.1462-2920.2000.00071.x 11243256 Parulekar N. N. Kolekar P. Jenkins A. Kleiven S. Utkilen H. Johansen A. (2017). Characterization of bacterial community associated with phytoplankton bloom in a eutrophic lake in South Norway using 16S rRNA gene amplicon sequence analysis. PLoS One 12:e0173408. 10.1371/journal.pone.0173408 28282404 Porter T. M. Hajibabaei M. (2018). Scaling up: a guide to high-throughput genomic approaches for biodiversity analysis. Mol. Ecol. 27 313338. 10.1111/mec.14478 29292539 Quiblier C. Wood S. Echenique-Subiabre I. Heath M. Villeneuve A. Humbert J.-F. (2013). A review of current knowledge on toxic benthic freshwater cyanobacteria - Ecology, toxin production and risk management. Water Res. 47 54645479. 10.1016/j.watres.2013.06.042 23891539 R Core Team, (2016). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. R Studio Team, (2015). RStudio: Integrated Development for R. Boston, MA: RStudio, Inc. Seifert M. McGregor G. Eaglesham G. Wickramasinghe W. Shaw G. (2007). First evidence for the production of cylindrospermopsin and deoxy-cylindrospermopsin by the freshwater benthic cyanobacterium, Lyngbya wollei (Farlow ex Gomont) Speziale and Dyck. Harmf. Algae 6 7380. 10.1016/j.hal.2006.07.001 Su X. Steinman A. D. Tang X. Xue Q. Zhao Y. Xie L. (2017). Response of bacterial communities to cyanobacterial harmful algal blooms in Lake Taihu, China. Harmf. Algae 68 168177. 10.1016/j.hal.2017.08.007 28962977 Wickham H. (2016). ggplot2: Elegant Graphics for Data Analysis. New York, NY: Springer-Verlag. Wood S. Puddick J. (2017). The abundance of toxic genotypes is a key contributor to anatoxin variability in Phormidium-dominated benthic mats. Mar. Drugs 15:307. 10.3390/md15100307 29019928 Wood S. A. Atalah J. Wagenhoff A. Brown L. Doehring K. Young R. G. (2017a). Effect of river flow, temperature, and water chemistry on proliferations of the benthic anatoxin-producing cyanobacterium Phormidium. Freshw. Sci. 36 6376. 10.1086/690114 Wood S. A. Depree C. Brown L. McAllister T. Hawes I. (2015). Entrapped sediments as a source of phosphorus in epilithic cyanobacterial proliferations in low nutrient rivers. PLoS One 10:e0141063. 10.1371/journal.pone.0141063 26479491 Wood S. A. Heath M. W. Kuhajek J. Ryan K. G. (2010). Fine-scale spatial variability in anatoxin-a and homoanatoxin-a concentrations in benthic cyanobacterial mats: implication for monitoring and management. J. Appl. Microbiol. 109 20112018. 10.1111/j.1365-2672.2010.04831.x 20722874 Wood S. A. Selwood A. I. Rueckert A. Holland P. T. Milne J. R. Smith K. F. (2007). First report of homoanatoxin-a and associated dog neurotoxicosis in New Zealand. Toxicon 50 292301. 10.1016/j.toxicon.2007.03.025 17517427 Wood S. A. Smith F. M. J. Heath M. W. Palfroy T. Gaw S. Young R. G. (2012). Within-mat variability in anatoxin-a and homoanatoxin-a production among benthic Phormidium (Cyanobacteria) strains. Toxins 4 900912. 10.3390/toxins4100900 23162704 Wood S. A. Puddick J. Fleming R. Heussner A. H. (2017b). Detection of anatoxin-producing Phormidium in a New Zealand farm pond and an associated dog death. New Zealand J. Bot. 55 3646. 10.1080/0028825x.2016.1231122 Woodhouse J. N. Kinsela A. S. Collins R. N. Bowling L. C. Honeyman G. L. Holliday J. K. (2015). Microbial communities reflect temporal changes in cyanobacterial composition in a shallow ephemeral freshwater lake. ISME J. 10:1337. 10.1038/ismej.2015.218 26636552 Yilmaz P. Parfrey L. W. Yarza P. Gerken J. Pruesse E. Quast C. (2014). The SILVA and “All-species Living Tree Project (LTP)” taxonomic frameworks. Nucleic Acids Res. 42 D643D648. 10.1093/nar/gkt1209 24293649 Zeglin L. H. (2015). Stream microbial diversity in response to environmental changes: review and synthesis of existing research. Front. Microbiol. 6:454. 10.3389/fmicb.2015.00454 26042102 Zinger L. Bonin A. Alsos I. G. Bálint M. Bik H. Boyer F. (2019). DNA metabarcoding—Need for robust experimental designs to draw sound ecological conclusions. Mol. Ecol. 28 18571862. 10.1111/mec.15060 31033079
      ‘Oh, my dear Thomas, you haven’t heard the terrible news then?’ she said. ‘I thought you would be sure to have seen it placarded somewhere. Alice went straight to her room, and I haven’t seen her since, though I repeatedly knocked at the door, which she has locked on the inside, and I’m sure it’s most unnatural of her not to let her own mother comfort her. It all happened in a moment: I have always said those great motor-cars shouldn’t be allowed to career about the streets, especially when they are all paved with cobbles as they are at Easton Haven, which are{331} so slippery when it’s wet. He slipped, and it went over him in a moment.’ My thanks were few and awkward, for there still hung to the missive a basting thread, and it was as warm as a nestling bird. I bent low--everybody was emotional in those days--kissed the fragrant thing, thrust it into my bosom, and blushed worse than Camille. "What, the Corner House victim? Is that really a fact?" "My dear child, I don't look upon it in that light at all. The child gave our picturesque friend a certain distinction--'My husband is dead, and this is my only child,' and all that sort of thing. It pays in society." leave them on the steps of a foundling asylum in order to insure [See larger version] Interoffice guff says you're planning definite moves on your own, J. O., and against some opposition. Is the Colonel so poor or so grasping—or what? Albert could not speak, for he felt as if his brains and teeth were rattling about inside his head. The rest of[Pg 188] the family hunched together by the door, the boys gaping idiotically, the girls in tears. "Now you're married." The host was called in, and unlocked a drawer in which they were deposited. The galleyman, with visible reluctance, arrayed himself in the garments, and he was observed to shudder more than once during the investiture of the dead man's apparel. HoME香京julia种子在线播放 ENTER NUMBET 0016lkchain.com.cn
      www.gnsuzu.com.cn
      hhccgo.com.cn
      giiild.com.cn
      www.iholdzhu.net.cn
      edjruu.com.cn
      uqsboc.com.cn
      www.pdwjhu.com.cn
      www.tinuan.com.cn
      www.odchsl.com.cn
      处女被大鸡巴操 强奸乱伦小说图片 俄罗斯美女爱爱图 调教强奸学生 亚洲女的穴 夜来香图片大全 美女性强奸电影 手机版色中阁 男性人体艺术素描图 16p成人 欧美性爱360 电影区 亚洲电影 欧美电影 经典三级 偷拍自拍 动漫电影 乱伦电影 变态另类 全部电 类似狠狠鲁的网站 黑吊操白逼图片 韩国黄片种子下载 操逼逼逼逼逼 人妻 小说 p 偷拍10幼女自慰 极品淫水很多 黄色做i爱 日本女人人体电影快播看 大福国小 我爱肏屄美女 mmcrwcom 欧美多人性交图片 肥臀乱伦老头舔阴帝 d09a4343000019c5 西欧人体艺术b xxoo激情短片 未成年人的 插泰国人夭图片 第770弾み1 24p 日本美女性 交动态 eee色播 yantasythunder 操无毛少女屄 亚洲图片你懂的女人 鸡巴插姨娘 特级黄 色大片播 左耳影音先锋 冢本友希全集 日本人体艺术绿色 我爱被舔逼 内射 幼 美阴图 喷水妹子高潮迭起 和后妈 操逼 美女吞鸡巴 鸭个自慰 中国女裸名单 操逼肥臀出水换妻 色站裸体义术 中国行上的漏毛美女叫什么 亚洲妹性交图 欧美美女人裸体人艺照 成人色妹妹直播 WWW_JXCT_COM r日本女人性淫乱 大胆人艺体艺图片 女同接吻av 碰碰哥免费自拍打炮 艳舞写真duppid1 88电影街拍视频 日本自拍做爱qvod 实拍美女性爱组图 少女高清av 浙江真实乱伦迅雷 台湾luanlunxiaoshuo 洛克王国宠物排行榜 皇瑟电影yy频道大全 红孩儿连连看 阴毛摄影 大胆美女写真人体艺术摄影 和风骚三个媳妇在家做爱 性爱办公室高清 18p2p木耳 大波撸影音 大鸡巴插嫩穴小说 一剧不超两个黑人 阿姨诱惑我快播 幼香阁千叶县小学生 少女妇女被狗强奸 曰人体妹妹 十二岁性感幼女 超级乱伦qvod 97爱蜜桃ccc336 日本淫妇阴液 av海量资源999 凤凰影视成仁 辰溪四中艳照门照片 先锋模特裸体展示影片 成人片免费看 自拍百度云 肥白老妇女 女爱人体图片 妈妈一女穴 星野美夏 日本少女dachidu 妹子私处人体图片 yinmindahuitang 舔无毛逼影片快播 田莹疑的裸体照片 三级电影影音先锋02222 妻子被外国老头操 观月雏乃泥鳅 韩国成人偷拍自拍图片 强奸5一9岁幼女小说 汤姆影院av图片 妹妹人艺体图 美女大驱 和女友做爱图片自拍p 绫川まどか在线先锋 那么嫩的逼很少见了 小女孩做爱 处女好逼连连看图图 性感美女在家做爱 近距离抽插骚逼逼 黑屌肏金毛屄 日韩av美少女 看喝尿尿小姐日逼色色色网图片 欧美肛交新视频 美女吃逼逼 av30线上免费 伊人在线三级经典 新视觉影院t6090影院 最新淫色电影网址 天龙影院远古手机版 搞老太影院 插进美女的大屁股里 私人影院加盟费用 www258dd 求一部电影里面有一个二猛哥 深肛交 日本萌妹子人体艺术写真图片 插入屄眼 美女的木奶 中文字幕黄色网址影视先锋 九号女神裸 和骚人妻偷情 和潘晓婷做爱 国模大尺度蜜桃 欧美大逼50p 西西人体成人 李宗瑞继母做爱原图物处理 nianhuawang 男鸡巴的视屏 � 97免费色伦电影 好色网成人 大姨子先锋 淫荡巨乳美女教师妈妈 性nuexiaoshuo WWW36YYYCOM 长春继续给力进屋就操小女儿套干破内射对白淫荡 农夫激情社区 日韩无码bt 欧美美女手掰嫩穴图片 日本援交偷拍自拍 入侵者日本在线播放 亚洲白虎偷拍自拍 常州高见泽日屄 寂寞少妇自卫视频 人体露逼图片 多毛外国老太 变态乱轮手机在线 淫荡妈妈和儿子操逼 伦理片大奶少女 看片神器最新登入地址sqvheqi345com账号群 麻美学姐无头 圣诞老人射小妞和强奸小妞动话片 亚洲AV女老师 先锋影音欧美成人资源 33344iucoom zV天堂电影网 宾馆美女打炮视频 色五月丁香五月magnet 嫂子淫乱小说 张歆艺的老公 吃奶男人视频在线播放 欧美色图男女乱伦 avtt2014ccvom 性插色欲香影院 青青草撸死你青青草 99热久久第一时间 激情套图卡通动漫 幼女裸聊做爱口交 日本女人被强奸乱伦 草榴社区快播 2kkk正在播放兽骑 啊不要人家小穴都湿了 www猎奇影视 A片www245vvcomwwwchnrwhmhzcn 搜索宜春院av wwwsee78co 逼奶鸡巴插 好吊日AV在线视频19gancom 熟女伦乱图片小说 日本免费av无码片在线开苞 鲁大妈撸到爆 裸聊官网 德国熟女xxx 新不夜城论坛首页手机 女虐男网址 男女做爱视频华为网盘 激情午夜天亚洲色图 内裤哥mangent 吉沢明歩制服丝袜WWWHHH710COM 屌逼在线试看 人体艺体阿娇艳照 推荐一个可以免费看片的网站如果被QQ拦截请复制链接在其它浏览器打开xxxyyy5comintr2a2cb551573a2b2e 欧美360精品粉红鲍鱼 教师调教第一页 聚美屋精品图 中韩淫乱群交 俄罗斯撸撸片 把鸡巴插进小姨子的阴道 干干AV成人网 aolasoohpnbcn www84ytom 高清大量潮喷www27dyycom 宝贝开心成人 freefronvideos人母 嫩穴成人网gggg29com 逼着舅妈给我口交肛交彩漫画 欧美色色aV88wwwgangguanscom 老太太操逼自拍视频 777亚洲手机在线播放 有没有夫妻3p小说 色列漫画淫女 午间色站导航 欧美成人处女色大图 童颜巨乳亚洲综合 桃色性欲草 色眯眯射逼 无码中文字幕塞外青楼这是一个 狂日美女老师人妻 爱碰网官网 亚洲图片雅蠛蝶 快播35怎么搜片 2000XXXX电影 新谷露性家庭影院 深深候dvd播放 幼齿用英语怎么说 不雅伦理无需播放器 国外淫荡图片 国外网站幼幼嫩网址 成年人就去色色视频快播 我鲁日日鲁老老老我爱 caoshaonvbi 人体艺术avav 性感性色导航 韩国黄色哥来嫖网站 成人网站美逼 淫荡熟妇自拍 欧美色惰图片 北京空姐透明照 狼堡免费av视频 www776eom 亚洲无码av欧美天堂网男人天堂 欧美激情爆操 a片kk266co 色尼姑成人极速在线视频 国语家庭系列 蒋雯雯 越南伦理 色CC伦理影院手机版 99jbbcom 大鸡巴舅妈 国产偷拍自拍淫荡对话视频 少妇春梦射精 开心激动网 自拍偷牌成人 色桃隐 撸狗网性交视频 淫荡的三位老师 伦理电影wwwqiuxia6commqiuxia6com 怡春院分站 丝袜超短裙露脸迅雷下载 色制服电影院 97超碰好吊色男人 yy6080理论在线宅男日韩福利大全 大嫂丝袜 500人群交手机在线 5sav 偷拍熟女吧 口述我和妹妹的欲望 50p电脑版 wwwavtttcon 3p3com 伦理无码片在线看 欧美成人电影图片岛国性爱伦理电影 先锋影音AV成人欧美 我爱好色 淫电影网 WWW19MMCOM 玛丽罗斯3d同人动画h在线看 动漫女孩裸体 超级丝袜美腿乱伦 1919gogo欣赏 大色逼淫色 www就是撸 激情文学网好骚 A级黄片免费 xedd5com 国内的b是黑的 快播美国成年人片黄 av高跟丝袜视频 上原保奈美巨乳女教师在线观看 校园春色都市激情fefegancom 偷窥自拍XXOO 搜索看马操美女 人本女优视频 日日吧淫淫 人妻巨乳影院 美国女子性爱学校 大肥屁股重口味 啪啪啪啊啊啊不要 操碰 japanfreevideoshome国产 亚州淫荡老熟女人体 伦奸毛片免费在线看 天天影视se 樱桃做爱视频 亚卅av在线视频 x奸小说下载 亚洲色图图片在线 217av天堂网 东方在线撸撸-百度 幼幼丝袜集 灰姑娘的姐姐 青青草在线视频观看对华 86papa路con 亚洲1AV 综合图片2区亚洲 美国美女大逼电影 010插插av成人网站 www色comwww821kxwcom 播乐子成人网免费视频在线观看 大炮撸在线影院 ,www4KkKcom 野花鲁最近30部 wwwCC213wapwww2233ww2download 三客优最新地址 母亲让儿子爽的无码视频 全国黄色片子 欧美色图美国十次 超碰在线直播 性感妖娆操 亚洲肉感熟女色图 a片A毛片管看视频 8vaa褋芯屑 333kk 川岛和津实视频 在线母子乱伦对白 妹妹肥逼五月 亚洲美女自拍 老婆在我面前小说 韩国空姐堪比情趣内衣 干小姐综合 淫妻色五月 添骚穴 WM62COM 23456影视播放器 成人午夜剧场 尼姑福利网 AV区亚洲AV欧美AV512qucomwwwc5508com 经典欧美骚妇 震动棒露出 日韩丝袜美臀巨乳在线 av无限吧看 就去干少妇 色艺无间正面是哪集 校园春色我和老师做爱 漫画夜色 天海丽白色吊带 黄色淫荡性虐小说 午夜高清播放器 文20岁女性荫道口图片 热国产热无码热有码 2015小明发布看看算你色 百度云播影视 美女肏屄屄乱轮小说 家族舔阴AV影片 邪恶在线av有码 父女之交 关于处女破处的三级片 极品护士91在线 欧美虐待女人视频的网站 享受老太太的丝袜 aaazhibuo 8dfvodcom成人 真实自拍足交 群交男女猛插逼 妓女爱爱动态 lin35com是什么网站 abp159 亚洲色图偷拍自拍乱伦熟女抠逼自慰 朝国三级篇 淫三国幻想 免费的av小电影网站 日本阿v视频免费按摩师 av750c0m 黄色片操一下 巨乳少女车震在线观看 操逼 免费 囗述情感一乱伦岳母和女婿 WWW_FAMITSU_COM 偷拍中国少妇在公车被操视频 花也真衣论理电影 大鸡鸡插p洞 新片欧美十八岁美少 进击的巨人神thunderftp 西方美女15p 深圳哪里易找到老女人玩视频 在线成人有声小说 365rrr 女尿图片 我和淫荡的小姨做爱 � 做爱技术体照 淫妇性爱 大学生私拍b 第四射狠狠射小说 色中色成人av社区 和小姨子乱伦肛交 wwwppp62com 俄罗斯巨乳人体艺术 骚逼阿娇 汤芳人体图片大胆 大胆人体艺术bb私处 性感大胸骚货 哪个网站幼女的片多 日本美女本子把 色 五月天 婷婷 快播 美女 美穴艺术 色百合电影导航 大鸡巴用力 孙悟空操美少女战士 狠狠撸美女手掰穴图片 古代女子与兽类交 沙耶香套图 激情成人网区 暴风影音av播放 动漫女孩怎么插第3个 mmmpp44 黑木麻衣无码ed2k 淫荡学姐少妇 乱伦操少女屄 高中性爱故事 骚妹妹爱爱图网 韩国模特剪长发 大鸡巴把我逼日了 中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片 大胆女人下体艺术图片 789sss 影音先锋在线国内情侣野外性事自拍普通话对白 群撸图库 闪现君打阿乐 ady 小说 插入表妹嫩穴小说 推荐成人资源 网络播放器 成人台 149大胆人体艺术 大屌图片 骚美女成人av 春暖花开春色性吧 女亭婷五月 我上了同桌的姐姐 恋夜秀场主播自慰视频 yzppp 屄茎 操屄女图 美女鲍鱼大特写 淫乱的日本人妻山口玲子 偷拍射精图 性感美女人体艺木图片 种马小说完本 免费电影院 骑士福利导航导航网站 骚老婆足交 国产性爱一级电影 欧美免费成人花花性都 欧美大肥妞性爱视频 家庭乱伦网站快播 偷拍自拍国产毛片 金发美女也用大吊来开包 缔D杏那 yentiyishu人体艺术ytys WWWUUKKMCOM 女人露奶 � 苍井空露逼 老荡妇高跟丝袜足交 偷偷和女友的朋友做爱迅雷 做爱七十二尺 朱丹人体合成 麻腾由纪妃 帅哥撸播种子图 鸡巴插逼动态图片 羙国十次啦中文 WWW137AVCOM 神斗片欧美版华语 有气质女人人休艺术 由美老师放屁电影 欧美女人肉肏图片 白虎种子快播 国产自拍90后女孩 美女在床上疯狂嫩b 饭岛爱最后之作 幼幼强奸摸奶 色97成人动漫 两性性爱打鸡巴插逼 新视觉影院4080青苹果影院 嗯好爽插死我了 阴口艺术照 李宗瑞电影qvod38 爆操舅母 亚洲色图七七影院 被大鸡巴操菊花 怡红院肿么了 成人极品影院删除 欧美性爱大图色图强奸乱 欧美女子与狗随便性交 苍井空的bt种子无码 熟女乱伦长篇小说 大色虫 兽交幼女影音先锋播放 44aad be0ca93900121f9b 先锋天耗ばさ无码 欧毛毛女三级黄色片图 干女人黑木耳照 日本美女少妇嫩逼人体艺术 sesechangchang 色屄屄网 久久撸app下载 色图色噜 美女鸡巴大奶 好吊日在线视频在线观看 透明丝袜脚偷拍自拍 中山怡红院菜单 wcwwwcom下载 骑嫂子 亚洲大色妣 成人故事365ahnet 丝袜家庭教mp4 幼交肛交 妹妹撸撸大妈 日本毛爽 caoprom超碰在email 关于中国古代偷窥的黄片 第一会所老熟女下载 wwwhuangsecome 狼人干综合新地址HD播放 变态儿子强奸乱伦图 强奸电影名字 2wwwer37com 日本毛片基地一亚洲AVmzddcxcn 暗黑圣经仙桃影院 37tpcocn 持月真由xfplay 好吊日在线视频三级网 我爱背入李丽珍 电影师傅床戏在线观看 96插妹妹sexsex88com 豪放家庭在线播放 桃花宝典极夜著豆瓜网 安卓系统播放神器 美美网丝袜诱惑 人人干全免费视频xulawyercn av无插件一本道 全国色五月 操逼电影小说网 good在线wwwyuyuelvcom www18avmmd 撸波波影视无插件 伊人幼女成人电影 会看射的图片 小明插看看 全裸美女扒开粉嫩b 国人自拍性交网站 萝莉白丝足交本子 七草ちとせ巨乳视频 摇摇晃晃的成人电影 兰桂坊成社人区小说www68kqcom 舔阴论坛 久撸客一撸客色国内外成人激情在线 明星门 欧美大胆嫩肉穴爽大片 www牛逼插 性吧星云 少妇性奴的屁眼 人体艺术大胆mscbaidu1imgcn 最新久久色色成人版 l女同在线 小泽玛利亚高潮图片搜索 女性裸b图 肛交bt种子 最热门有声小说 人间添春色 春色猜谜字 樱井莉亚钢管舞视频 小泽玛利亚直美6p 能用的h网 还能看的h网 bl动漫h网 开心五月激 东京热401 男色女色第四色酒色网 怎么下载黄色小说 黄色小说小栽 和谐图城 乐乐影院 色哥导航 特色导航 依依社区 爱窝窝在线 色狼谷成人 91porn 包要你射电影 色色3A丝袜 丝袜妹妹淫网 爱色导航(荐) 好男人激情影院 坏哥哥 第七色 色久久 人格分裂 急先锋 撸撸射中文网 第一会所综合社区 91影院老师机 东方成人激情 怼莪影院吹潮 老鸭窝伊人无码不卡无码一本道 av女柳晶电影 91天生爱风流作品 深爱激情小说私房婷婷网 擼奶av 567pao 里番3d一家人野外 上原在线电影 水岛津实透明丝袜 1314酒色 网旧网俺也去 0855影院 在线无码私人影院 搜索 国产自拍 神马dy888午夜伦理达达兔 农民工黄晓婷 日韩裸体黑丝御姐 屈臣氏的燕窝面膜怎么样つぼみ晶エリーの早漏チ○ポ强化合宿 老熟女人性视频 影音先锋 三上悠亚ol 妹妹影院福利片 hhhhhhhhsxo 午夜天堂热的国产 强奸剧场 全裸香蕉视频无码 亚欧伦理视频 秋霞为什么给封了 日本在线视频空天使 日韩成人aⅴ在线 日本日屌日屄导航视频 在线福利视频 日本推油无码av magnet 在线免费视频 樱井梨吮东 日本一本道在线无码DVD 日本性感诱惑美女做爱阴道流水视频 日本一级av 汤姆avtom在线视频 台湾佬中文娱乐线20 阿v播播下载 橙色影院 奴隶少女护士cg视频 汤姆在线影院无码 偷拍宾馆 业面紧急生级访问 色和尚有线 厕所偷拍一族 av女l 公交色狼优酷视频 裸体视频AV 人与兽肉肉网 董美香ol 花井美纱链接 magnet 西瓜影音 亚洲 自拍 日韩女优欧美激情偷拍自拍 亚洲成年人免费视频 荷兰免费成人电影 深喉呕吐XXⅩX 操石榴在线视频 天天色成人免费视频 314hu四虎 涩久免费视频在线观看 成人电影迅雷下载 能看见整个奶子的香蕉影院 水菜丽百度影音 gwaz079百度云 噜死你们资源站 主播走光视频合集迅雷下载 thumbzilla jappen 精品Av 古川伊织star598在线 假面女皇vip在线视频播放 国产自拍迷情校园 啪啪啪公寓漫画 日本阿AV 黄色手机电影 欧美在线Av影院 华裔电击女神91在线 亚洲欧美专区 1日本1000部免费视频 开放90后 波多野结衣 东方 影院av 页面升级紧急访问每天正常更新 4438Xchengeren 老炮色 a k福利电影 色欲影视色天天视频 高老庄aV 259LUXU-683 magnet 手机在线电影 国产区 欧美激情人人操网 国产 偷拍 直播 日韩 国内外激情在线视频网给 站长统计一本道人妻 光棍影院被封 紫竹铃取汁 ftp 狂插空姐嫩 xfplay 丈夫面前 穿靴子伪街 XXOO视频在线免费 大香蕉道久在线播放 电棒漏电嗨过头 充气娃能看下毛和洞吗 夫妻牲交 福利云点墦 yukun瑟妃 疯狂交换女友 国产自拍26页 腐女资源 百度云 日本DVD高清无码视频 偷拍,自拍AV伦理电影 A片小视频福利站。 大奶肥婆自拍偷拍图片 交配伊甸园 超碰在线视频自拍偷拍国产 小热巴91大神 rctd 045 类似于A片 超美大奶大学生美女直播被男友操 男友问 你的衣服怎么脱掉的 亚洲女与黑人群交视频一 在线黄涩 木内美保步兵番号 鸡巴插入欧美美女的b舒服 激情在线国产自拍日韩欧美 国语福利小视频在线观看 作爱小视颍 潮喷合集丝袜无码mp4 做爱的无码高清视频 牛牛精品 伊aⅤ在线观看 savk12 哥哥搞在线播放 在线电一本道影 一级谍片 250pp亚洲情艺中心,88 欧美一本道九色在线一 wwwseavbacom色av吧 cos美女在线 欧美17,18ⅹⅹⅹ视频 自拍嫩逼 小电影在线观看网站 筱田优 贼 水电工 5358x视频 日本69式视频有码 b雪福利导航 韩国女主播19tvclub在线 操逼清晰视频 丝袜美女国产视频网址导航 水菜丽颜射房间 台湾妹中文娱乐网 风吟岛视频 口交 伦理 日本熟妇色五十路免费视频 A级片互舔 川村真矢Av在线观看 亚洲日韩av 色和尚国产自拍 sea8 mp4 aV天堂2018手机在线 免费版国产偷拍a在线播放 狠狠 婷婷 丁香 小视频福利在线观看平台 思妍白衣小仙女被邻居强上 萝莉自拍有水 4484新视觉 永久发布页 977成人影视在线观看 小清新影院在线观 小鸟酱后丝后入百度云 旋风魅影四级 香蕉影院小黄片免费看 性爱直播磁力链接 小骚逼第一色影院 性交流的视频 小雪小视频bd 小视频TV禁看视频 迷奸AV在线看 nba直播 任你在干线 汤姆影院在线视频国产 624u在线播放 成人 一级a做爰片就在线看狐狸视频 小香蕉AV视频 www182、com 腿模简小育 学生做爱视频 秘密搜查官 快播 成人福利网午夜 一级黄色夫妻录像片 直接看的gav久久播放器 国产自拍400首页 sm老爹影院 谁知道隔壁老王网址在线 综合网 123西瓜影音 米奇丁香 人人澡人人漠大学生 色久悠 夜色视频你今天寂寞了吗? 菲菲影视城美国 被抄的影院 变态另类 欧美 成人 国产偷拍自拍在线小说 不用下载安装就能看的吃男人鸡巴视频 插屄视频 大贯杏里播放 wwwhhh50 233若菜奈央 伦理片天海翼秘密搜查官 大香蕉在线万色屋视频 那种漫画小说你懂的 祥仔电影合集一区 那里可以看澳门皇冠酒店a片 色自啪 亚洲aV电影天堂 谷露影院ar toupaizaixian sexbj。com 毕业生 zaixian mianfei 朝桐光视频 成人短视频在线直接观看 陈美霖 沈阳音乐学院 导航女 www26yjjcom 1大尺度视频 开平虐女视频 菅野雪松协和影视在线视频 华人play在线视频bbb 鸡吧操屄视频 多啪啪免费视频 悠草影院 金兰策划网 (969) 橘佑金短视频 国内一极刺激自拍片 日本制服番号大全magnet 成人动漫母系 电脑怎么清理内存 黄色福利1000 dy88午夜 偷拍中学生洗澡磁力链接 花椒相机福利美女视频 站长推荐磁力下载 mp4 三洞轮流插视频 玉兔miki热舞视频 夜生活小视频 爆乳人妖小视频 国内网红主播自拍福利迅雷下载 不用app的裸裸体美女操逼视频 变态SM影片在线观看 草溜影院元气吧 - 百度 - 百度 波推全套视频 国产双飞集合ftp 日本在线AV网 笔国毛片 神马影院女主播是我的邻居 影音资源 激情乱伦电影 799pao 亚洲第一色第一影院 av视频大香蕉 老梁故事汇希斯莱杰 水中人体磁力链接 下载 大香蕉黄片免费看 济南谭崔 避开屏蔽的岛a片 草破福利 要看大鸡巴操小骚逼的人的视频 黑丝少妇影音先锋 欧美巨乳熟女磁力链接 美国黄网站色大全 伦蕉在线久播 极品女厕沟 激情五月bd韩国电影 混血美女自摸和男友激情啪啪自拍诱人呻吟福利视频 人人摸人人妻做人人看 44kknn 娸娸原网 伊人欧美 恋夜影院视频列表安卓青青 57k影院 如果电话亭 avi 插爆骚女精品自拍 青青草在线免费视频1769TV 令人惹火的邻家美眉 影音先锋 真人妹子被捅动态图 男人女人做完爱视频15 表姐合租两人共处一室晚上她竟爬上了我的床 性爱教学视频 北条麻妃bd在线播放版 国产老师和师生 magnet wwwcctv1024 女神自慰 ftp 女同性恋做激情视频 欧美大胆露阴视频 欧美无码影视 好女色在线观看 后入肥臀18p 百度影视屏福利 厕所超碰视频 强奸mp magnet 欧美妹aⅴ免费线上看 2016年妞干网视频 5手机在线福利 超在线最视频 800av:cOm magnet 欧美性爱免播放器在线播放 91大款肥汤的性感美乳90后邻家美眉趴着窗台后入啪啪 秋霞日本毛片网站 cheng ren 在线视频 上原亚衣肛门无码解禁影音先锋 美脚家庭教师在线播放 尤酷伦理片 熟女性生活视频在线观看 欧美av在线播放喷潮 194avav 凤凰AV成人 - 百度 kbb9999 AV片AV在线AV无码 爱爱视频高清免费观看 黄色男女操b视频 观看 18AV清纯视频在线播放平台 成人性爱视频久久操 女性真人生殖系统双性人视频 下身插入b射精视频 明星潜规测视频 mp4 免賛a片直播绪 国内 自己 偷拍 在线 国内真实偷拍 手机在线 国产主播户外勾在线 三桥杏奈高清无码迅雷下载 2五福电影院凸凹频频 男主拿鱼打女主,高宝宝 色哥午夜影院 川村まや痴汉 草溜影院费全过程免费 淫小弟影院在线视频 laohantuiche 啪啪啪喷潮XXOO视频 青娱乐成人国产 蓝沢润 一本道 亚洲青涩中文欧美 神马影院线理论 米娅卡莉法的av 在线福利65535 欧美粉色在线 欧美性受群交视频1在线播放 极品喷奶熟妇在线播放 变态另类无码福利影院92 天津小姐被偷拍 磁力下载 台湾三级电髟全部 丝袜美腿偷拍自拍 偷拍女生性行为图 妻子的乱伦 白虎少妇 肏婶骚屄 外国大妈会阴照片 美少女操屄图片 妹妹自慰11p 操老熟女的b 361美女人体 360电影院樱桃 爱色妹妹亚洲色图 性交卖淫姿势高清图片一级 欧美一黑对二白 大色网无毛一线天 射小妹网站 寂寞穴 西西人体模特苍井空 操的大白逼吧 骚穴让我操 拉好友干女朋友3p