Can metabarcoding resolve intraspecific genetic diversity changes to environmental stressors? A test case using river macrozoobenthos

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 straight-forward and usually limited to one or few individual 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 mitochondrial genetic 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 looking at haplotype diversity. Here, stressor responses were only 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 river systems, only OTUs could be used, which were present in all systems. As OTU composition differed strongly between the rivers, this led to the exclusion of a high number of OTUs, especially in diverse river systems of good quality, which potentially diminished the genetic diversity patterns. To better understand responses of genetic diversity to environmental stressors for example in river ecosystems, it would be important to increase OTU overlap between sites of comparisons, e.g. by sampling a narrower stressor gradient, and to perform calibrated studies controlling for the number and individual genotypes. However, this pioneer study shows that the extraction of haplotypes from DNA metabarcoding datasets is a promising tool to simultaneously assess mitochondrial genetic diversity changes in response to environmental impacts for a metacommunity.


Introduction 67
Degradation, pollution, and exploitation of freshwater ecosystems has resulted in a drastic decline of 68 biodiversity (Vörösmarty et al., 2010;WWF, 2018). The magnitude of biodiversity loss depends on 69 stressor intensities as well as on resistance and resilience of the biotic communities (Dobson et al., 2006; 70 Elmqvist et al., 2003;Vörösmarty et al., 2010). So far, degradation and recovery processes have mostly 71 been studied at the level of species diversity (alpha diversity). However, the underlying genetic diversity 72 within species is an essential variable to consider in this context, as it determines the evolutionary 73 capacity of a species to adapt to changing environments. Genetic diversity in ecosystems is influenced 74 by four major processes: mutation, drift, migration, and selection (Vellend and Geber, 2005). A high  97 Sturmbauer et al., 1999), the level of intraspecific genetic diversity still goes unnoticed, when distance-98 based clustering is used to analyse the sequencing data. As an alternative, bioinformatic denoising 99 approaches can be used to obtain 'exact sequence variants' (ESVs) (

104
(COI) is used for DNA barcoding and a shorter part of this for metabarcoding (Hebert et al., 2003).

105
While the use of haploid, maternally-inherited mitochondrial markers has limitations for detailed 106 population genetic analyses (Ballard and Whitlock, 2004;Leese and Held, 2011), its utility to infer 107 insights into geographic structure and population diversity, and thereby ecological processes acting at 108 local or regional level, has often been demonstrated (e.g. Pauls    conditions and stressor impacts are also present in the Ennepe, and include near-natural rural sites, urban 127 sites stressed through occasional stormwater retention basin overflow, and sites with sewage treatment 128 plant inflow. In comparison, the river Sieg is considered as a stable, near-natural river system with a 129 good ecological and chemical status. Punctual stressor inflow is present through rainwater retention 130 basins but not immediately at sampling sites. By comparing communities between the different streams,

131
we want to test, if species and genetic diversity is correlated with present stressor gradients. Following 132 predictions from the ecological habitat concept, which links presence and abundance of species over 133 time to the available resources (see e.g. Van Dyck, 2012), we expect OTU diversity to be highest in the 134 river Sieg, moderate in river Ennepe, and lowest in river Emscher. According to the genetic erosion 135 hypothesis, we predict haplotype diversity to covary with OTU richness. We expect highest haplotype 136 number and diversity at the river Sieg due to its long-time stable good ecological conditions supporting 137 large and stable population sizes. Lower values are expected at the river Ennepe, where communities 138 are regularly affected by organic stressor influx and thus recurrent population decline, resulting in higher 139 genetic drift. The lowest values are expected for river Emscher due to the complete erasure of MZB diversity in the history of this river system caused by the usage as sewage transport system, and the still 141 prevalent stress level in many parts. We predict different patterns for taxa sensitive (EPT -142 Ephemeroptera, Plecoptera, Trichoptera) or resistant to organic pollution (PR -'Pollution Resistant' -143 Arhynchobdellida, Enchytraeida, Haplotaxida, Isopoda, Rhynchobdellida). Specifically, we assume that 144 OTU and haplotype diversity for EPT taxa will decline with increasing stress because of declining 145 population sizes, and eventual local species extinction. In contrast we expect that PR taxa will show 146 opposing patterns, because of their resistance to organic pollution and their ability to rather use organic

214
Haplotypes with similarity < 95 % to a deposited sequence in the database were excluded from further 215 analysis to prevent incorrect assignments potentially leading to the assessment of erroneous diversity 216 patterns. Read numbers per haplotype of technical PCR replicates were fused and the average was  e.g. E5) for EPT taxa, and 8 -257 haplotypes clustered into 4 -55 OTUs for PR taxa. As no plecopterans 238 were found at river Emscher, only ET taxa could be analysed for this river. Total OTU and haplotype 239 number of EPT taxa was affected by the river system with more OTUs and haplotypes at Ennepe and 240 Sieg than at Emscher (p < 0.001). Sample sites E4 and E5 showed remarkably high OTU and haplotype 241 numbers assigned to PR taxa compared to all other sample sites. However, no effect of the river system 242 was detected on PR taxa (haplotype number p = 0.09) (Fig. 2). Number of counted individuals before 243 laboratory processing differed between all river systems (p < 0.001) and seasons (p < 0.05) (Table S2).

244
Within streams, individual numbers did not significantly differ between sampling sites. However, no 245 correlation was detected between total specimen number per site and season, and average haplotype 246 number per OTU (Fig. S1). shared OTU when comparing all three river systems (EPT: p < 0.05, PR: p < 0.05) (Fig. 3) (Fig. 3).

288
Further, the comparison of shared OTUs between the three river systems revealed an effect of river 289 system on nucleotide diversity for EPT taxa (p < 0.05). Average nucleotide diversity of taxa was higher 290 at river Ennepe (0.00215) and Sieg (0.00266) than at the Emscher (0.00104). In contrast, PR taxa showed 291 a higher nucleotide diversity at the Emscher (0.00143) than at the Ennepe (0.00215), but only when 292 comparing OTUs shared between those two rivers (p < 0.05, Fig. S6

299
We plotted total OTU number assigned to EPT and PR taxa against the average haplotype number per 300 OTU and sample site for the four datasets of shared OTUs (Fig. 4 A-H)

303
Emscher to Sieg and Ennepe (Fig. 4 A). A weaker correlation was observed for EPT taxa comparing 304 Ennepe and Sieg (Fig. 4 D). In addition, correlations were significant for PR taxa comparing Emscher,

308
Emscher and Ennepe (Fig. 4 B) and Emscher and Sieg (Fig. 4 C). The separation was most distinct 309 between Emscher and Ennepe. Separation of river systems due to total number of PR taxa (x axis Fig.   310 E-H) was less distinct than separation based on total ET taxa.

Discussion 314
We expected a general effect of stressors on alpha (OTU) diversity and intraspecific genetic diversity in 315 the three river systems. In line with these expectations, the heavily impacted river Emscher showed the 316 lowest OTU number compared with the other two systems. However, no significant differences between 317 Ennepe and Sieg were detected, and an even higher number of highly pollution-sensitive stoneflies 318 (Plecoptera) was found at the stronger impacted river Ennepe. Stressor effects at this river therefore 319 seem to be less intense than expected with no obvious effects on MZB communities. As expected, 320 haplotype numbers per sample site were lower in river Emscher than in Emscher and Sieg, especially Emscher and Sieg and correlations between total number of PR taxa and genetic diversity.

344
The data presented hold great potential to inform on metacommunity and metapopulation structure and 345 processes. However, several limitations need to be considered when further testing and applying the

386
Individual sample size: Lastly, while the strength of the approach is that thousands of specimens were 387 processed at once, there is little control about the individual number of specimens per species or OTU.

388
Comparisons of genetic diversity rely on the number of specimens sampled and future studies should 389 carefully control for this in order to test the reliability and robustness of the approach.

Conclusion 392
Our study shows that metabarcoding-derived haplotype information can be utilized to assess effects of 393 environmental variables on intraspecific genetic diversity using the example of three differently-

394
impacted German river systems. Even if the choice of data filtering thresholds is still a trade-off between 395 rare 'real' sequences and erroneous, artificially generated ones, we were able to link stressor effects to 396 genetic variability of aquatic macroinvertebrate communities. Sites with good and stable ecological 397 conditions showed higher genetic diversity than stressed sites, which is coupled with OTU diversity.

398
However, due to a low OTU overlap between river systems, genetic diversity analyses were based only 399 on subsets, including all shared OTUs. This subsampling induced the exclusion of variability and 400 ecological specialists, especially at highly diverse sample sites and might skew actual differences.