Environmental DNA time series analysis of a temperate stream reveals distinct seasonal community and functional shifts

Environmental DNA (eDNA) extracted from water is routinely used in river biodiversity research, and via metabarcoding eDNA can provide comprehensive taxa lists with little effort and cost. However, eDNA‐based species detection in streams and rivers may be influenced by sampling season and other key factors such as water temperature and discharge. Research linking these factors and also informing on the potential of eDNA metabarcoding to detect shifts in ecological signatures, such as species phenology and functional feeding groups across seasons, is missing. To address this gap, we collected water samples every 2 weeks for 15 months at a long‐term ecological research (LTER) site and at three different positions in the river's cross section, specifically the water surface, riverbed, and riverbank. For these 102 samples, we analyzed macroinvertebrate species and molecular operational taxonomic unit (OTU) richness and temporal community turnover across seasons based on cytochrome c oxidase subunit I (COI) metabarcoding data. Using Generalized Additive Models, we found a significant influence of sampling season on species richness. Community turnover followed a cyclic pattern, reflecting the continuous change of the macroinvertebrate community throughout the year (“seasonal clock”). Although water temperature had no influence on the inferred species richness, higher discharge reduced the number of Annelida and Ephemeroptera species detectable with eDNA. Most macroinvertebrate taxa showed the highest species richness in spring, in particular merolimnic species with univoltine life cycles. Further, we detected an increase in the proportion of shredders in winter and parasites in summer. Our results show the usefulness of highly resolved eDNA metabarcoding time series data for ecological research and biodiversity monitoring in streams and rivers.

Biodiversity in streams and rivers is being impacted by multiple anthropogenic stressors (Jackson et al., 2016).To understand these impacts, their functional consequences, and management effectiveness taxonomically highly resolved information with high temporal resolution is important.However, such information is difficult to obtain through traditional morphological assessments as many invertebrate species are small or present only in juvenile stages that are difficult to identify.Molecular taxonomic approaches, in particular metabarcoding of environmental DNA (eDNA) collected from water, offer a fast and cost-effective way to assess biodiversity and are routinely used in aquatic bioassessments around the world (Deiner et al., 2017).Environmental DNA metabarcoding is based on extracted DNA shed by organisms via sloughed cells, feces, gametes, or other particles into the water and is thus a non-invasive method to assess community composition because assessment is based on water rather than organismal bulk samples (Taberlet, Coissac, Hajibabaei et al., 2012;Taberlet, Coissac, Pompanon et al., 2012;Valentini et al., 2016).DNA metabarcoding uses high-throughput sequencing methods to generate comprehensive taxa lists (Brantschen et al., 2022;Leese et al., 2021).However, since the reference databases used to assign taxonomic names to the obtained sequences are still incomplete (Weigand et al., 2019), not all sequences can be assigned to species level.Therefore, molecular operational taxonomic units (OTUs) that are generated according to genetic distance-based similarity thresholds can be used as surrogates for species.Using OTUs in addition to species can reveal further insights into ecological processes (e.g., Beermann et al., 2018).
Despite the obvious advantages, several factors hinder the direct interpretation of eDNA data (Barnes & Turner, 2016;Harrison et al., 2019).First, several abiotic factors can influence DNA transport and detectability and may thus distort the inferred community (Barnes & Turner, 2016;Harrison et al., 2019), such as discharge and water temperature.Discharge is an important factor influencing eDNA detectability in streams because high discharge could lead to more species being detected by eDNA signals from transported DNA or whole organisms being swept downstream (Carraro et al., 2018;Fremier et al., 2019;Shogren et al., 2017).In contrast, high discharge can also dilute the eDNA signal thus making it more difficult to detect all present species (e.g., Thalinger et al., 2021), which may particularly impede the detection of rare species that are already at low abundance.In addition to discharge, temperature also affects eDNA detectability (Strickler et al., 2015), either negatively if higher temperatures reduce eDNA persistence due to increased enzymatic activity or positively if higher temperatures increase DNA shedding rates (Jo et al., 2019;Kasai et al., 2020;Strickler et al., 2015).These potentially contrasting effects of discharge and temperature make it difficult to predict how they will affect estimates of community composition.
Second, as a consequence of the phenology of organisms, eDNA detectability in streams might also be influenced by sampling season.
Depending on the stream type and region, characteristic abundance patterns can be found for different macroinvertebrate orders, genera, and species throughout the year and across years (Cowell et al., 2004;Wagner, 2004;Wagner et al., 2011).In addition, the biology of the different macroinvertebrate taxa has a strong effect on seasonal community composition.One important factor is differences in organismal life cycles.While hololimnic species (species with a fully aquatic life cycle) are presumed to be present in the water the whole year, merolimnic species (species with aquatic larvae and aerial adults) leave the water after hatching and have distinct emergence periods, often lasting up to few months, which can lead to a sudden decline in sampled benthic communities (Füreder et al., 2005;Jackson & Füreder, 2006).
Besides life cycle-based community composition changes, streams also differ systematically with respect to their functional feeding groups (FFGs) in both space and time (Vannote et al., 1980).For example, shredders are typically more abundant in autumn, when the amount of allochthonous material in streams is highest (Cummins et al., 1989) and grazers in spring and summer, due to sun exposure supporting the growth of large biofilms.The functional composition of macroinvertebrate communities in the form of different FFGs affects ecosystem functioning and is therefore also included in bioassessment approaches ( Šporka et al., 2006).Detecting these community dynamics patterns is important in aquatic ecology.There is ample evidence that these types of seasonal differences are reflected in eDNA metabarcoding data (Bista et al., 2017;Dunn et al., 2017;Zizka et al., 2020) and that season or even month of sampling lead to different biological assessment results with eDNA metabarcoding (Jensen et al., 2021;Zizka et al., 2020).
While biodiversity studies addressing larger temporal scales often suffer from an insufficient resolution (Jackson & Fuereder, 2006;Pilotto et al., 2020), studies using high-resolution temporal data at smaller scales are still scarce.Especially for macroinvertebrates, eDNA metabarcoding data has the potential to assess small temporal changes in community composition time-and costeffectively to complement future long-term bioassessment of streams.
Using time series data from a long-term ecological research (LTER; Mirtl et al., 2018) site, the aim of this study was to test the effect of discharge and temperature, and sampling season on stream macroinvertebrate community composition and species richness determined from eDNA.Further, we included three biological replicates to test for small-scale differences between sampling positions in the river's cross section.Small-scale differences in sampling position may recover different lotic communities (Macher & Leese, 2017) but alternatively, sampling position may have no effect given that eDNA can be transported over long distances of >12 km in streams (Deiner & Altermatt, 2014), which may homogenize eDNA community signals across sampling positions (Macher et al., 2021).
The time series comprises in total 102 eDNA samples taken every 2 weeks for 15 months (from 24 May 2017 to 29 August 2018) at three sampling positions at the same location which served as replicates: (i) the river surface; (ii) the river bottom; and (iii) the river bank.
Community composition and species richness were determined using high-throughput mitochondrial cytochrome c oxidase subunit I (COI) gene metabarcoding.We tested three hypotheses: 1. Changes in species richness and community composition will be driven by seasonality and differences in seasons will follow a cyclic pattern throughout the year.
2. Differences in species richness will reflect the diverging life cycles of different mero-and hololimnic taxa.For example, the biomass of merolimnic taxa like Ephemeroptera, Plecoptera, and Trichoptera (EPT) will decline in summer after emergence leading to lower species richness detectable with eDNA metabarcoding, while differences will be less pronounced for hololimnic taxa like Annelida and Coleoptera.
3. eDNA metabarcoding will detect seasonal differences in community composition for different FFGs.In particular, more grazer species will be present in spring and summer, whereas more shredder species will be found in autumn and winter based on seasonal food availability.For parasite species that are dependent on other organisms, overall species occurrence will not follow any seasonal pattern, as it is linked to the presence of different host taxa.
Moreover, we also tested the effects of discharge and water temperature on species richness, although we had no a priori expectations of what these effects would be.

| Sampling, filtration, and extraction
Sampling, filtration, and extraction methods are described in detail in Leese et al. (2021).In brief, we sampled in total 102 water samples from the 86 km long Kinzig River every 2 weeks within a LTER (Mirtl et al., 2018) site at the Rhine-Main-Observatory (RMO, https://deims.org/9f9ba137-342d-4813-ae58-a60911c3abc1) from 24 May 2017 to 29 August 2018.We took three 1-L samples per sampling day at the surface in the middle of the stream, 10 cm above the stream bed in the middle of the stream, with varying stream depths between 0.48 and 2.56 m, and at the riverbank (Figure 1).After sampling, water samples were put on ice in the field and later stored in a freezer at À20 C until further processing.Water samples were filtered using a vacuum pump and DNA captured on 0.45 μm cellulose nitrate membrane filters (Nalgene, diameter 47 mm).For each day and person who filtered, a negative filter control was included, where only the surrounding air was filtered to control for air and filter contamination.
Filters were stored in 1.5 mL Eppendorf tubes filled with 96% denatured ethanol.DNA extraction was performed following the salt-EtOH-precipitation protocol reported in Weiss and Leese (2016).Subsequent steps included RNA digestion using 1.5 μL RNase A (10 mg/ mL, Thermo Fisher Scientific) per sample and clean-up of the samples using a Qiagen MinElute Kit to obtain high-quality DNA.

| DNA metabarcoding
A two-step PCR approach with two replicates per extract was applied (Zizka et al., 2019).For the first step, four length-varying primers for the primer pair fwhF2/EPTDr2n (Leese et al., 2021;Vamos et al., 2017) were used, with each of them having a universal tail attached.Length variation was due to inline shifts (0-3 Ns) between the universal tail and the primer sequence to maximize diversity between sequences (Elbrecht & Leese, 2015).For the second step, primers matching the universal tail with an i5/i7 index and P5/P7 Illumina adapter attached were used (Buchner, Beermann et al., 2021;Buchner, Haase et al., 2021).We applied an improved PCR protocol compared to the previous study (Leese et al., 2021) The library was sequenced on a single flow cell on an Illumina HiSeq X sequencer using the 300 cycle kit (150 bp paired-end reads) with 5% Phi-X spike-in added by the sequencing company Macrogen Europe B.V. Ten of the 102 samples of the dataset were amplified and sequenced separately in another library, as they had already been processed in another study (Leese et al., 2021).For bioinformatics analysis, sequences from the 10 previously sequenced samples (Short-Read Archive, Project number: PRJNA664693) and the sequences from the remaining 92 samples (deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB61051) were processed together.

| Bioinformatics
Sequence reads were processed as stated in our previous study in detail (Leese, et al., 2021).Briefly, JAMP v0.67 (https://github.com/VascoElbrecht/JAMP; Elbrecht et al., 2018) was used on default settings to merge paired-end reads and, where needed, to build the reverse complements of the sequences.Primer sequences were removed.To retain only reads of the expected fragment length, sequences with a deviation of >15 bp were excluded from further analyses.Reads with an expected maximum error of >0.5 and singletons were removed before clustering the sequences with a similarity ≥97% to OTUs.To maximize the number of reads retained, the dereplicated sequences, including singletons, were mapped with a similarity of ≥97% to the generated OTU dataset.Only OTUs with a minimal read abundance of 0.01% in at least one sample were retained for further analyses.OTU centroid sequences were compared to the BOLD database for taxonomic annotation using BOL-Digger 1.1.4(Buchner & Leese, 2020).For further analyses, we only considered OTUs with a similarity of ≥90% to a reference sequence in BOLD.OTUs with a similarity of ≥98% were assigned to species, ≥95% to genus, and ≥90% to family level.Replicates were merged with reads summed up and divided by two for each OTU.OTUs for which conflicting taxonomic results were found were checked manually, taking into account if reference specimens were identified by taxonomic experts.Further, the obtained taxa list was compared to the RMO database, which contains detailed information on morphologically identified taxa occurring in this area, and the taxa list was additionally checked by three taxonomic experts to exclude terrestrial taxa and taxa that are impossible or unlikely to occur in the study area.

| Species trait data
We split our dataset into merolimnic and hololimnic species because we expected to see differences in species and OTU richness over our sampling period based on their different life cycles.Additionally, we looked more in detail into species and OTU richness of the most abundant merolimnic (EPT and Diptera) and hololimnic orders (Annelida and Coleoptera) and their associated FFGs.To assign life cycle characteristics and FFGs to our taxa, we used trait data from freshwaterecology.info (Schmidt-Kloiber & Hering, 2015) for a given taxon if available.If not available, we used data from a higher taxonomic level, such as genus-level trait values if species-level values were missing.For taxa with no trait values at the species-, genus-, or family-level, we used an averaging procedure whereby values were averaged across species within the genus if species-level values were available and across all taxa within the family if not (Kunz et al., 2022).We used the proportion of the different FFGs as the proportion of the overall detected feeding types divided by the number of all species found per sampling day.

| Environmental data
Data for discharge was accessible from a gaging station close to Hanau (www.hlnug.de)<1 km downstream from the sampling site with no tributaries in between.Discharge was measured daily with minimum, maximum, and average discharge data available.Water temperature data were recorded hourly over the entire sampling period at the sampling site using a temperature logger.

| Statistical analysis
All statistical analyses were performed in R v4.0.5 Venn diagrams were constructed using http://bioinformatics.psb.ugent.be/webtools/Venn/.The effect of sampling day, sampling position, and temperature and discharge on overall species/OTU detection was tested using generalized additive models (GAM) with a Poisson distribution and a "log" link function implemented by the gam() function in the R package mgcv (Wood, 2015).Where overdispersion was detected, we corrected the standard errors using a quasi-GAM model.Sampling day was defined as a linear description of time starting with the first day of sampling and discharge and temperature were respectively quantified as the average discharge and mean hourly temperature of each sampling day (Table S1).To test for the effect of our predictor variable sampling position on our response variable OTU richness/ read abundance a Kruskal-Wallis rank sum test was performed, as our data were non-normal and non-linear and we had three group variables.Further, we divided our dataset into the two groups merolimnic and hololimnic taxa and subsetted these groups into the most abundant taxa, Diptera and EPT for merolimnic and Annelida and Coleoptera for hololimnic.Then, we used GAMs to determine the effect of sampling day, temperature, and discharge on species/OTU richness separately for each of the eight groups.To assess the effect of sampling day on the proportion for each FFGs we used a GAM with a binomial distribution and a "logit" link function.Only the term for sampling day was smoothed in all GAMS via cubic regression splines with a basis dimension k = 9.Model diagnostics were checked via plots of residuals versus fitted values (Figures S20-S31).For all models, the chosen k was checked using gam.check() and no model had a significant p-value.The marginal effects, the slopes of the prediction equation, of the significant explanatory variables were plotted at the mean of the covariates to incorporate the possible minimal effects that the insignificant explanatory variables could potentially have on the slope of the significant explanatory variables.
An ANOSIM test was used to test for the effect of sampling position on OTU and species temporal beta-diversity (Jaccard distance) defined as the shift in the identities of named taxa in a specified assemblage over two or more time points (Magurran et al., 2019).
Additionally, we partitioned temporal beta-diversity into its components turnover-reflecting the replacement of species between samples-and nestedness-occurring when species loss or gain causes species-poor samples to resemble a strict subset of species-rich samples (Baselga, 2010) using the betapart package (Baselga et al., 2018) and visualized their change over the time series as pairwise comparisons in relation to the first timepoint.For comparisons between seasons, we grouped our dataset into spring (March, April, and May), summer (June, July, and August), autumn (September, October, and November), and winter (December, January, and February) samples.
To visualize differences between sampling positions and seasons for species and OTU community composition, non-metric multidimensional scaling (NMDS) based on Jaccard's distance was conducted using the package vegan (Oksanen et al., 2022).
To identify species that are most strongly associated with a season or only found at a specific season, an indicator species analysis based on the presence-absence was performed using the R package indicspecies with the function multipatt() (De Cáceres et al., 2010).
The association of species to seasons was calculated using the indicator value index from Dufrêne & Legendre, 1997.We used the indicator value index instead of a correlation-based index to get results that best match the observed presence of a species rather than the preferences of one season over another.Using preferences could lead to the problem of species being reported as indicators for a specific season, although the species was also almost frequently present in other seasons.For biomonitoring aspects, it is therefore favorable to use the observed presence of species.The indicator value index ranges between 0 and 1 (with 0 indicating no association) and is calculated as the product of two quantities A, the probability that the surveyed site belongs to the target group (i.e., season) when the species has been found, and B, how frequently the species is found in the different samples from each season.The function determines which season corresponds to the highest association value for each species and then tests the statistical significance of this relationship using a permutation test.We corrected for size differences between groups, as we had different sample sizes between seasons, using the implemented value within the multipatt() function.

| RESULTS
We obtained 372,431,352 reads and 4702 OTUs for the 102 eDNA samples after quality filtering and excluding OTUs with read abundance <0.01%per sample.Reads among samples were unevenly distributed ranging from 88,446 to 4,451,399 with a mean read number of 1,880,966.We excluded two samples that had only 3282 and 7928 reads, respectively, to avoid unrepresentative biodiversity data due to too low sequencing depth.After the database and taxonomist checks, 2840 OTUs remained which were assigned to 506 species.The uneven read distribution was no concern as neither OTU richness nor species richness were correlated with read numbers (see Figures S1   and S2).Eight controls, including filter blanks, with no DNA were also included to test for contamination during field or laboratory proceedings.Their read numbers ranged from zero to four, indicating no substantial contamination.The vast majority of reads were assigned to metazoans, especially arthropods and annelids, with <1% of reads assigned to non-metazoan taxa.In all samples, chironomids had the highest read numbers and constituted 55% to 95% of every sample.

| Effect of sampling day, sampling position, and environment
There was a significant relationship between sampling day and both species (F = 5.933, p-value = <0.001,ndf (numerator degrees of freedom) = 4.055, ddf (denominator degrees of freedom) = 90.945,R 2 (adj.)= 0.336), and OTU (F = 2.765, p-value = 0.030, ndf = 3.32, ddf = 91.68,R 2 (adj.)= 0.103) richness.In contrast, neither temperature, discharge, nor sampling position had an effect on the number of species or OTUs detected.The highest species and OTU richness was detected in spring with the lowest richness in the summer months for all three sampling positions (Figures 2 and S3).The vast majority of macroinvertebrate species were detected at all three sampling positions (280), also with over 99% of the reads assigned to these (Figure 2).Although some species were exclusively found at one of the three sampling positions, the relative occurrence of these species was low in contrast to shared species (Figure S4).We found no significant differences between surface, riverbed, or riverbank samples based on OTU richness (Kruskal-Wallis rank sum test: χ 2 = 0.085, pvalue = 0.957, ndf = 2, ddf = 97).
Additionally, the NMDS analysis showed no differences in species or OTU beta-diversity between sampling positions based on ANOSIM tests ( p-value: 0.7371/0.7079)(Figures S5 and S6).
Based on these results, samples from the different sampling positions were merged for each sampling day and subsequent analyses focused on sampling time and seasonal differences.NMDS analysis for the merged samples revealed a strong impact of sampling season on community structure (Figures 3 and S7).
Samples from the following years cluster closely together, for example, 07 June 2017/06 June 2018 and 23 May 2018/24 May 2017.Summer and winter samples were separated on the first axis, spring and autumn samples on the second axis, producing a cyclic temporal beta-diversity pattern resembling a "seasonal clock."When temporal beta-diversity was partitioned into turnover and nestedness, turnover was the driving factor for change in beta-diversity while nestedness was of less importance (Figure S8).Over the year, turnover increased over the winter samples and decreased again in the summer with nestedness showing the opposite trend, reaching its highest value in summer 2018 (Figure S8).
Temperature had no effect on the detection of species nor OTU richness for any of the tested groups.Higher discharge negatively affected the detection of Ephemeroptera and overall hololimnic and Annelida species (E: χ 2 = 4.213, p-value = 0.04, ndf = 1, ddf = 28.292,R 2 (adj.)= 0.384; hololimnic: χ 2 = 6.639, pvalue = 0.01, ndf = 1, ddf = 25.502,R 2 (adj.)= 0.439; Annelida: S11 and S12).Per 10 m 3 /s increase in discharge Ephemeroptera lost 1 species, hololimnic taxa 5 species and Annelida were congruent with the NMDS and occurrence plots showing a clear distinction between summer and winter samples.In some cases, more than one OTU was assigned to the same species.Here, most OTUs showed similar association values.However, for Baetis vernus the different OTUs showed different associations even though it had no significant association with any season(s) when all OTUs were combined.
One OTU was significantly associated with autumn, one with autumn + spring and two with autumn + spring + summer.

| Influence of sampling time on species richness and community composition
Our results indicated that differences in species richness were caused by seasonality but not sample replicates, which matched our expectations.High similarity among sampling positions likely occurred due to the complete mixing of eDNA in the water body with smaller differences occurring between replicates, the three water samples taken within the sampling site, probably due to sampling error, and not by systematic differences between positions (Macher et al., 2021).The stream sampled in our study had a water depth <1 m and sampling positions were <10 m apart with a mostly turbulent flow, all of which can reduce community heterogeneity (see Fremier et al., 2019;Shogren et al., 2017 for further discussions).
Temporal beta diversity shifted in a seasonally cyclical and gradual pattern, which was primarily driven by turnover, with turnover and nestedness both returning to a similar level after 1 year.The potential of eDNA in detecting seasonal shifts in community composition has been already shown for marine environments (Jensen et al., 2022;Liu et al., 2022) but not for freshwater.

| Influence of different life cycles on species and OTU detection
We hypothesized that species richness will differ depending on life cycle differences between mero-and hololimnic taxa.Accordingly, we found a strong effect of sampling season on all tested groups, except Coleoptera.Species richness of merolimnic and hololimnic taxa were both affected by sampling day but in a different way.In contrast to merolimnic taxa, hololimnic taxa showed a more stable number of species between autumn 2017 and spring 2018 likely due to the absence of emergence events, as hypothesized.For the most species-rich hololimnic group, Annelida, most species were detected in autumn, likely due to an increase of annelid populations and hence stronger eDNA signal in the water.The highest densities of aquatic annelid species are known for summer and autumn (Cook, 1969;Learner et al., 1978;Smith, 1986) with densities for the species Dero digita peaking in autumn (Smith, 1986) and in our data set the species was only detected once in October.Contrary to our second hypothesis, we did detect a lower species richness in summer also for the hololimnic group Annelida, but the decrease was more steady instead of an abrupt drop in species richness after peaking as for merolimnic species.
In general, merolimnic species detection in the water depends on several factors, for example, time and duration of the flight period and emergence, the abundance of the species in the water, number of generations per year, larval development, and occurrence of dormancy/diapause.Diptera were the most species-rich merolimnic order with most taxa belonging to the family Chironomidae.We detected a high diversity of Diptera in spring and early summer which is consistent with the results of other studies, implying a reduced diversity in winter and the highest diversity in early summer before emergence, as the emergence of chironomids is correlated with high temperatures (Armitage et al., 2012;Bista et al., 2017).The high number of species in spring is likely due to a high growth rate preparing the larvae for emergence in summer, also explaining the drop in detected species richness for the later summer month.For the species-rich mero-and hololimnic groups, except Annelida, no effect of sampling day on OTU richness was detected, possibly due to the high number of OTUs per group which can show distinct responses to seasonal changes even within the same species.In contrast, for Ephemeroptera, Plecoptera, and Trichoptera, only OTU richness was affected by sampling day, likely due to the lower numbers of species and therefore smaller shifts in species richness over time.Plecoptera and Trichoptera OTU richness showed one peak each, in winter and early spring, which is likely linked to an increase in biomass that will probably reach its peak right before emergence coinciding with our second hypothesis expecting a low richness in summer.Plecoptera OTU richness peaked a bit earlier than Trichoptera OTU richness which is congruent with an on average earlier emergence of many Plecoptera species compared to most Trichoptera species (Graf et al., 2008;Graf et al., 2009;Graf, Lorenz et al., 2022, Graf, Murphy et al., 2022;Schmidt-Kloiber & Hering, 2015).By comparison, Ephemeroptera OTU richness had two peaks, one in autumn and one in late spring.Species details revealed that from 21 Ephemeroptera species, ten have a bivoltine (mainly genus Baetis), one a flexible (Caenis beskidensis), and two a semivoltine (genus Ephemera) lifecycle (Buffagni et al., 2009;Buffagni et al., 2022;Schmidt-Kloiber & Hering, 2015).In contrast, the detected Plecoptera and Trichoptera were mostly univoltine species.
Occurrence even between species within a merolimnic order differed between seasons.The indicator species analysis revealed that G. pellucidus is detectable in all winter and spring months and absent in summer which coincides with the species having a long flight period after emergence and a known dormancy (Graf et al., 2008;Graf, Lorenz et al., 2022, Graf, Murphy et al., 2022;Schmidt-Kloiber & Hering, 2015).In contrast, Philopotamus ludificatus has also a long flight period and dormancy but was detected less frequently but had similar read numbers as G. pellucidus (Graf et al., 2008;Graf, Lorenz et al., 2022, Graf, Murphy et al., 2022;Schmidt-Kloiber & Hering, 2015).As neither differences between life cycle characteristics nor read numbers were present, the less frequent detection of  (Buffagni et al., 2009;Buffagni et al., 2022;Schmidt-Kloiber & Hering, 2015) and therefore, larvae of B. rhodani are probably more frequent.Additionally, the OTUs detected belonging to species B. vernus differed in their seasonal occurrence with one OTU only occurring in autumn.It is known that different OTUs can show distinct responses to environmental changes (Beermann et al., 2018), and for the B. vernus group cryptic diversity has been recorded (Ståhls & Savolainen, 2008).This strengthens the assumption that the differences in detection between the OTUs are based on different responses to environmental changes and therefore an underestimation of the diversity within B. vernus.A limitation in using eDNA to assess patterns of seasonality is the persistence of eDNA in the environment for up to several days or weeks, and apart from that, the uncertainty of capturing living or dead cells (Dejean et al., 2011;Pilliod et al., 2014;Thomsen, Kielgast, Iversen, Møller, et al., 2012;Thomsen, Kielgast, Iversen, Wiuf et al., 2012).Nonetheless, our results demonstrate that the patterns we found are consistent with the phenology of the different taxa thus further encouraging that the DNA we detected mostly originated from living organisms.

| Influence of season on the distribution of feeding types
For two of the FFGs, we observed an influence of sampling time on the proportion of FFGs, partly confirming our third hypothesis.Shredders showed a large increase from summer to winter and a decrease afterward.As shredders feed on allochthonous material such as decomposed leaf litter, the increase in the proportion of shredders in winter is likely due to an increase in their food after fallen leaves accumulate in the stream (Cummins, 1974;Cummins et al., 1989).
Although we expected to detect a higher number of grazers in the warmer months, we did not detect an effect of sampling day on the proportion of grazers.As the number of detected species was in general extremely low in the summer samples, the lack of seasonal effect on the proportion of grazers detected is likely the result of the extremely dry and warm summer in 2018 reducing the abundance of temperature-sensitive taxa.Since we only looked at presence/ absence data, it is possible that abundance data would have shown higher proportions of grazers in summer.Proportion of parasites was highest in summer and their occurrence was likely linked to their hosts.The exclusively parasitic larvae of Sisyra terminalis, which was frequently found in the spring and summer but not winter samples, is known to parasitize on sponges and the species Spongilla lacustris and Eunapius fragilis, which were also present in our dataset, are known to form colonies during warm seasons and produce gemmules over the cold seasons (Gugel, 2001;Weißmair, 1994).Additionally, also the Trichoptera Ceraclea nigronervosa, which feeds on Spongillidae (Graf et al., 2002), was detected in the summer and spring months but not in the winter months, further encouraging the influence of the abundance of sponges on the detection of their hosts/predators.These results are promising for future eDNA-based biomonitoring, as we could detect shifts in the proportion of different FFGs and therefore changes in the functional composition of the macroinvertebrate community over time, which could help to determine changes in ecosystem functioning (Minshall et al., 1992).

| Influence of water temperature and discharge on species and OTU richness
We neither detected an influence of temperature on overall species nor on OTU richness, thus indicating that other factors like sampling time or season had a greater effect on eDNA species and OTU recovery.In contrast, we found a negative effect of increasing discharge on hololimnic taxa, Annelida, and Ephemeroptera species richness.For Ephemeroptera with a p-value close to 0.05 this should be interpreted as a trend and future studies are needed to verify the effect of discharge.For hololimnic species, where Annelida were the most species-rich group, a possible explanation could be that many Annelida specimens being swept away with high discharge is low as most annelids are embedded within the sediment and are potentially moving even deeper into the sediment with high discharge.This could lead to a higher dilution effect and therefore lower annelid DNA concentrations as specimens are not washed away.
In the summer months, an extremely low species and OTU richness was detected for all groups with August 2018 being the month with the lowest richness.The summer 2018 was outstandingly warm with low rainfall and high sun exposure.The exceptionally warm temperatures coupled with low rainfall probably affected the abundance and growth rates of the more sensitive species therefore causing a decrease in detection rates but also increased DNAse activity or microbial break DNA down leading to faster DNA degradation rates could be an influencing factor (Barnes & Turner, 2016;Nielsen et al., 2007).

| CONCLUSION
Our data show how eDNA-based time-series data on a stream macroinvertebrate community can support ecological research and biomonitoring of rivers, as we detected seasonal differences and small temporal shifts in community composition and key functional feeding groups.As a clear advantage of eDNA metabarcoding, our study detected many species of morphologically neglected taxa such as annelids and chironomids that are underrepresented in other studies.
By including a broader range of taxa in the analysis, eDNA metabarcoding could function as a complementary tool for long-term monitoring when assessing differences in macroinvertebrate community composition on a highly resolved temporal scale.
, PCR reactions for both PCR steps contained: 1Â Multiplex PCR Master Mix (Qiagen Multiplex PCR Plus Kit), 0.1 μM of each primer, 0.5 Â CoralLoad Dye, 1 μL DNA extract/unpurified PCR product (1st step/2nd step) filled up with Rnase-free water to a total volume of 25 μL.Cycling conditions were as follows: first step PCRs consisted of 5 min initial F I G U R E 1 eDNA sampling location within the long-term ecological research (LTER) site Rhine-main-observatory (RMO) in Hesse, Germany.Right: sampling location with the three sampling positions marked with a red cross.(a) sampling position: riverbank in winter, (b) zoomed-in sampling position: riverbank in summer, (c) sampling position: surface and riverbed in summer, and (d) zoomed-in sampling position: surface and riverbed in winter.[Color figure can be viewed at wileyonlinelibrary.com] denaturation at 95 C, followed by a touchdown of 10 cycles of 30 s at 95 C, 90 s at 64 C (À1 C each cycle), and 30 s at 72 C and additional 25 cycles of 30 s at 95 C, 90 s at 54 C, and 30 s at 72 C, followed by a final elongation for 10 min at 68 C. Second step PCRs consisted of 5 min initial denaturation at 95 C followed by 18 cycles of 30 s at 95 C and 2 min at 72 C, followed by a final elongation for 10 min at 68 C. Eight PCR negative controls were added to check for contamination during PCR and tag switching.PCR products were quantified and equimolarly pooled into one library before purification and left-sided size selection with a ratio of 0.76Â (Beckman Coulter).

F
I G U R E 2 Comparison of the three sampling positions.Left: Venn diagram for the different combinations of sampling positions.For each combination read proportion in relation to overall macroinvertebrate reads is given in percent, numbers of species shown in brackets.Right: species richness over time for each sampling position.Points represent pure species richness, solid lines represent the fitted values of the generalized additive models (GAM).[Color figure can be viewed at wileyonlinelibrary.com]F I G U R E 3 Non-metric multidimensional scaling (NMDS) based on Jaccard's distance for MZB species for the 34 samples after merging sampling position samples for each time point.Autumn samples: September, October, and November; Winter samples: December, January, and February; Spring samples: March, April, and May; Summer samples: June, July, and August.[Color figure can be viewed at wileyonlinelibrary.com] (115), and the lowest richness was detected in summer (60) (Figure S13), as well as for the most species-rich merolimnic subgroup Diptera (spring: 95 species and summer: 50 species) (Figure S14).While the detection of merolimnic subgroups Ephemeroptera, Plecoptera, and Trichoptera species richness was unaffected by sampling day, it had an effect on the detection of OTU (E: χ 2 = 39.11,p-value = <0.001,ndf = 7.27, ddf = 23.73,R 2 (adj.)= 0.433; P: χ 2 = 21.01,p-value = 0.005, ndf = 6.156, ddf = 24.844,R 2 (adj.)= 0.781; T: χ 2 = 20.11,p-value = 0.002, ndf = 4.636, ddf = 26.364,R 2 (adj.)= 0.738).For Ephemeroptera, most OTUs were detected in autumn (15) and in late spring (20).In contrast, the detected Plecoptera and Trichoptera OTU richness was highest from December to March (10) and January to March ( analysis identified most species indicative for a certain season or seasons when either winter or spring samples were included in the combinations: autumn + spring + winter: 35 indicator species; spring + winter: 30; winter: 24; spring: 20.Groups including the summer samples showed the lowest number of indicator species: autumn + summer: 2 indicator species; summer: 1; autumn + summer + winter: 1; spring + summer + winter: 1, indicating the presence of few species characteristic only for summer.Some species were exclusively found in one season, like the Plecoptera Protonemura nitida in autumn and the Trichoptera Chaetopteryx major in winter (see F I G U R E 4 Occurrence plot of the most common merolimnic species that occur at least in 10% of the samples.Diptera are excluded because species richness was too high to visualize clearly.A separate plot for Diptera is presented in Figure S19.[Color figure can be viewed at wileyonlinelibrary.com] also Figures 4 and S18).Species not found in summer, but quite frequently in all other seasons, were Trichoptera Silo piceus, Potamophylax latipennis, Plecoptera Nemoura flexuosa, Diptera Tipula paludosa, Rheocricotopus atripes, Coleoptera Platambus maculatus, and Annelida Tubifex tubifex.Species that were almost exclusively (high A value) and often (high B value) found in spring and winter included Trichoptera Glyphotaelius pellucidus, Limnephilus flavicornis, Halesus radiatus, and Plecoptera Nemurella pictetii, Isoperla grammatica.Here, results

F
I G U R E 5 Proportion of species belonging to a distinct functional feeding group over time.(a) shr (shredders), n = 120; (b) par (parasites), n = 7. Points represent pure values, solid lines represent the fitted values of the generalized additive models (GAM).[Color figure can be viewed at wileyonlinelibrary.com]

P
. ludificatus was likely because the species was rarer at the sampling location or because of other unknown factors influencing the detectability of the species.Other species like the Ephemeroptera Baetis rhodani and B. vernus are detectable throughout the year with B. vernus being less frequently detectable than B. rhodani.The more frequent detection of B. rhodani can be explained due to the species being also sometimes trivoltine and B. vernus being only known for a bivoltine life cycle