Research Article |
Corresponding author: Vera Marie Alida Zizka ( v.zizka@leibniz-zfmk.de ) Academic editor: Xin Zhou
© 2020 Vera Marie Alida Zizka, Martina Weiss, Florian Leese.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Zizka VMA, Weiss M, Leese F (2020) Can metabarcoding resolve intraspecific genetic diversity changes to environmental stressors? A test case using river macrozoobenthos. Metabarcoding and Metagenomics 4: e51925. https://doi.org/10.3897/mbmg.4.51925
|
Genetic diversity is the most basal level of biodiversity and determines the evolutionary capacity of species to adapt to changing environments, yet it is typically neglected in routine biomonitoring and stressor impact assessment. For a comprehensive analysis of stressor impacts on genetic diversity, it is necessary to assess genetic variants simultaneously in many individuals and species. Such an assessment is not as straightforward and usually limited to one or few focal species. However, nowadays species diversity can be assessed by analysing thousands of individuals of a community simultaneously with DNA metabarcoding. Recent bioinformatic advances also allow for the extraction of exact sequence variants (ESVs or haplotypes) in addition to Operational Taxonomic Units (OTUs). By using this new capability, we here evaluated if the analysis of intraspecific mitochondrial diversity in addition to species diversity can provide insights into responses of stream macrozoobenthic communities to environmental stressors. For this purpose, we analysed macroinvertebrate bulk samples of three German river systems with different stressor levels using DNA metabarcoding. While OTU and haplotype number were negatively correlated with stressor impact, this association was not as clear when studying haplotype diversity across all taxa. However, stressor responses were found for sensitive EPT (Ephemeroptera, Plecoptera, Trichoptera) taxa and those exceedingly resistant to organic stress. An increase in haplotype number per OTU and haplotype diversity of sensitive taxa was observed with an increase in ecosystem quality and stability, while the opposite pattern was detected for pollution resistant taxa. However, this pattern was less prominent than expected based on the strong differences in stressor intensity between sites. To compare genetic diversity among communities in river systems, we focussed on OTUs, which were present in all systems. As OTU composition differed strongly between rivers, this led to the exclusion of a high number of OTUs, especially in diverse river systems of good quality, which potentially diminished the increase in intraspecific diversity. To better understand responses of intraspecific genetic diversity to environmental stressors, for example in river ecosystems, it would be important to increase OTU overlap between compared sites, e.g. by sampling a narrower stressor gradient, and to perform calibrated studies controlling for the number of individuals and their haplotypes. However, this pioneer study shows that the extraction of haplotypes from DNA metabarcoding datasets is a promising source of information to simultaneously assess intraspecific diversity changes in response to environmental impacts for a metacommunity.
denoising, environmental impact, genetic erosion, haplotypes, pollution
Degradation, pollution, and exploitation of freshwater ecosystems have resulted in a drastic decline of biodiversity (
Nowadays, alpha diversity can be assessed with great resolution using DNA metabarcoding (
In our study, we wanted to further evaluate the potential of ESVs as a proxy for COI haplotypes in addition to OTU data obtained from environmental bulk sample metabarcoding, to analyse the impact of stressors on macrozoobenthos (MZB) communities in river ecosystems. MZB organisms play a key role in freshwater ecosystem functionality and include a wide range of taxonomic groups with often narrow and specific demands with respect to habitat conditions (
Macroinvertebrates were sampled according to Water Framework Directive compliant protocol (Meier et al. 2006) at six sites in the rivers Emscher and Sieg in autumn 2016 and 2017, and spring 2017 and 2018 (Figure
Samples were examined under a binocular (Leica S6E) to separate individuals from substrate. Substrate was discarded and individuals were counted and separated into two size categories (size class A: ≤ 25 mm, size class B: ≥ 25 mm) (see Elbrecht et al. 2017b for the procedure). Individuals of the two size classes were dried in petri dishes overnight and homogenised to fine powder with an IKA Ultra Turrax Tube Disperser (BMT-20-S-M sterile tubes, full speed for 30 min). Two times half of a spatula (~20 mg) was transferred into two Eppendorf tubes, 600 µl TNES buffer and 15 µl Proteinase K (10 mg/ml) were added per tube, and incubated overnight at 36 °C shaking at 250 rpm in a Thermoshaker (ThermoMixer C, Eppendorf). A salt extraction protocol (after Sunnucks and Hales 1996, adjusted as in
Sequences were analysed with JAMP-0.67 (https://github.com/VascoElbrecht/JAMP) including demultiplexing of data, paired-end-merging, and primer trimming, following standard settings. For haplotype extraction only reads of expected fragment length (421 bp) were included. A strict quality filtering was applied (maximal expected error max_ee = 0.3) and reads with an abundance < 0.003 % in a sample were excluded from the dataset. The algorithm Unoise3 (
After denoising and abundance filtering, on average 29,063 reads were present in Emscher samples, 36,289 reads in Sieg samples, and 51,981 reads in Ennepe samples. Because filtering thresholds were based on relative abundances (see Material and methods, Data analysis), reads per sample were not adjusted to uniform numbers. Samples contained 228–694 haplotypes, which clustered into 70–155 OTUs. OTU and haplotype number was higher at river Ennepe and Sieg than at river Emscher (p < 0.01, Fig.
Total OTU (upper part) and haplotype number (lower part) summed across all seasons for all river systems (Em 1-6 – Emscher, En 1-7 – Ennepe 1-6, S – Sieg). A) All benthic macroinvertebrate taxa assigned to a reference sequence in BOLD with > 95% sequence identity; B) EPT taxa (Ephemeroptera, Plecoptera, Trichoptera); C) PR taxa (Pollution Resistant : Enchytraeida, Haplotaxida, Isopoda, Rhynchobdellida).
To compare average haplotype number per OTU and haplotype diversity between river systems, we searched for OTUs present in most samples. The five most common OTUs, all occurring at more than 50 % of the analysed samples, were: OTU 2 (Baetis rhodani, 63 %), OTU 9 (Asellus aquaticus, 52 %), OTU 12 (Stylodrilus heringianus, 63 %), OTU 65 (Esolus parallelepipedus, 55 %), and OTU 107 (Microtendipes pedellus, 52 %). To further increase the number of shared OTUs between sites, samples collected at different seasons were merged (E1-E6, En1-En7, S1-S6). By this, we identified six OTUs occurring in more than 80 % of all sites (OTU 2: Baetis rhodani, 84 %; OTU 5: Orthocladius sp. A, 89 %, OTU 9: Asellus aquaticus, 84 %; OTU 45: Orthocladius sp. B, 84 %; OTU 67: Tanytarsus eminulus, 84 %, OTU 107: Microtendipes pedellus, 95 %). To increase the number of OTUs for analyses, all OTUs present in at least one of the samples per river system were included, resulting in four different datasets (Em-En-S: 78 shared OTUs, Em-En: 110 shared OTUs, Em-S: 125 shared OTUs, En-S: 155 shared OTUs). Per dataset > 47 % of shared OTUs were assigned to dipterans, of which the majority (> 90 %) were chironomids (Suppl. material
A) Average haplotype number per OTU for the four datasets of shared OTUs. B) Haplotype diversity for the four datasets of shared OTUs. Sensitive (EPT) and pollution resistant (‘PR’) taxa are shown. * indicates significant difference between regarded groups.
Further, the comparison of shared OTUs between the three river systems revealed an effect of river system on nucleotide diversity for EPT taxa (p < 0.05). Average nucleotide diversity of taxa was higher at river Ennepe (0.00215) and Sieg (0.00266) than at Emscher (0.00104). In contrast, PR taxa showed a higher nucleotide diversity at the Emscher (0.00143) than at the Ennepe (0.00215), but only when comparing OTUs shared between those two rivers (p < 0.05, Suppl. material
We plotted total OTU number assigned to EPT and PR taxa against the average haplotype number per OTU and sample site for the four datasets of shared OTUs (Fig.
We expected a general effect of stressors on alpha (OTU) diversity and intraspecific genetic diversity in the three river systems. In line with these expectations, the heavily impacted river Emscher showed the lowest OTU number compared with the other two systems. However, no significant differences between Ennepe and Sieg were detected, and an even higher number of highly pollution-sensitive stoneflies (Plecoptera) was found at the stronger impacted river Ennepe. Since strong stressor impact differences between river systems were expected, the lack of MZB community impacts between river Sieg and Ennepe seemed unexpected. However, further system or population specific factors, which were not specifically tested in the present study, could have influenced the observed pattern. These factors should be considered in the future (e.g. population history), also including individual numbers of investigated species. As expected, haplotype numbers per sample site were lower in river Emscher than in Ennepe and Sieg, especially for pollution sensitive mayfly and caddisfly taxa. Vice versa, pollution resistant (‘PR’) taxa had higher haplotype numbers in the river Emscher. Genetic diversity estimates inferred via average haplotype number per OTU, haplotype diversity as well as nucleotide diversity revealed higher values at river Sieg and Ennepe compared to the Emscher, but again no differences between the former two rivers and therefore supported results based on OTU and haplotype number. In this case, the hypothesis that sustainable good ecological conditions and stability at river Sieg induced stable and large population size, favouring a high level of genetic diversity in sensitive MZB communities, was supported (
The data presented, hold great potential to inform on metacommunity and metapopulation structure and processes. However, several limitations need to be considered when further testing and applying the approach:
OTU overlap: The fact, that no significant differences between average haplotype number per OTU and haplotype diversity were found for the other comparisons of E(P)T taxa (pairwise comparisons of all rivers), probably results from the small overlap of OTUs shared between rivers. The highly heterogeneous, site-specific stressor levels and variable time since restoration at Emscher further inflated variation in OTU composition, thereby limiting statistical power for comparisons of genetic variation based on highly frequent OTUs. We decided to reduce our dataset to OTUs at least present in one sample site of compared river systems to circumvent differences in intraspecific genetic variation due to species specific traits and sensitivities (e.g. Sturmbaue et al. 1999;
Mitochondrial single-gene marker: Beside the limitations due to a small overlap in OTUs between river systems, the utility of the mitochondrial COI marker as a measure of intraspecific genetic diversity and variability could have deflated signal strength. Even if various studies have successfully applied this marker for even more specific population analyses, the use of only a single gene can be misleading – and mitochondrial genes are especially unique due to their haploid structure, the lack of recombination and purely maternal inheritance (
Bioinformatic haplotype extraction: Furthermore, as outlined in previous studies (
Individual sample size: Lastly, while the strength of the approach is that thousands of specimens were processed at once, there is little control about the individual number of specimens per species or OTU. Comparisons of genetic diversity rely on the number of specimens sampled and future studies should carefully control for this in order to test the reliability and robustness of the approach.
Using macrozoobenthic taxa from three differently impacted German river systems, our study shows that denoised metabarcoding data can provide valuable information to assess effects of environmental variables on intraspecific genetic diversity. Even if the choice of data filtering thresholds is still a trade-off between rare ‘real’ sequences and erroneous, artificially generated ones, we were able to link stressor effects to changes of intraspecific genetic variation in aquatic macroinvertebrate communities. Sites with good and stable ecological conditions showed higher intraspecific diversity than stressed sites, which is also coupled with higher OTU diversity. However, due to a low OTU overlap between river systems, genetic diversity analyses were based only on subsets, including all shared OTUs. This subsampling induced the exclusion of variability and ecological specialists, especially at highly diverse sample sites and might have skewed actual differences. Due to these limitations in the underlying data, we cannot disentangle effects of stressors from e.g. population and colonisation and dynamics. Future studies need to address these aspects in more detail and should include several replicates of similar conditions, presuppose a threshold of overlapping OTUs in compared river systems and control for the number of individuals per OTU to allow for a statistically accurate comparison. However, in conclusion, our study adds to the growing number of studies that highlight the potential to extract haplotype information from metabarcoding datasets for more holistic biodiversity assessments.
This study was supported by the German Federal Ministry of Education and Research project “German Barcode of Life II” (GBOLII) to FL (FKZ 01LI1101 and 01LI1501). We want to thank the Emschergenossenschaft/Lippeverband and Ennepe-Ruhr-Kreis for support during this study and Bianca Peinert for her help in the field and the laboratory. FL and VMAZ are supported by the COST Action CA15219 DNAqua-Net.
Figure S1. Total number of aquatic macroinvertebrate individuals per sample and season plotted against the average haplotype number per OTU. Different colours indicate the three river systems
Data type: genetic
Figure S2. Four different datasets including shared OTUs between the different river systems (Emscher-Ennepe-Sieg, Emscher-Ennepe, Emscher-Sieg, Sieg-Ennepe). Number of OTUs is illustrated with taxonomic assignment on order level
Data type: genetic, occurance
Figure S3. Average haplotype number per OTU for the four different datasets of shared OTUs. Values are illustrated for all sample sites including all shared OTUs
Data type: genetic
Figure S4. Average haplotype number per OTU for the four different datasets of shared OTUs. Datasets are split into EPT (Ephemeroptera, Plecoptera, Trichoptera) and PR (‘Pollution Resistant’) taxa
Data type: genetic
Figure S5. Average haplotype diversity for all four datasets of shared OTUs seperated according to sample sites and EPT (Ephemeroptera, Plecoptera, Trichoptera) and PR (‘Pollution Resistant’) taxa
Data type: (measurement/occurence/multimedia/etc.)
Figure S6. Average nucleotide diversity for all four datasets of shared OTUs seperated according to sample sites and EPT (Ephemeroptera, Plecoptera, Trichoptera) and PR (‘Pollution Resistant’) taxa
Data type: genetic
Figure S7 – part 1. Haplotype network of the two most frequent EPT (Ephemeroptera, Plecoptera, Trichoptera) and PR (‘Pollution Resistant’) taxa
Data type: genetic
Table S1. Number of macroinvertebrate individuals per sample and season
Data type: individual counts