Environmental DNA provides greater insight to biodiversity and ecosystem function compared to traditional approaches, via spatio-temporal nestedness and turnover partitioning

Rapidly assessing biodiversity is essential for environmental monitoring, however traditional approaches are limited in the scope needed for most ecological systems. Environmental DNA (eDNA) based assessment offers increased scope, while reducing cost and time, compared to traditional methods. Here we investigated the effects of landuse and seasonality on headwater community richness and functional diversity, to assess spatio-temporal dynamics between eDNA and traditional methods. Environmental DNA provided greater biodiversity resolution with both methods resulting in complementary ndings. Community richness was seasonally linked, peaking in spring and summer, with temporal turnover having a greater effect on community composition compared to localized nestedness. Our assessment of ecosystem function shows community formation is driven by regional resource availability, implying regional management requirements. Our ndings show that eDNA based ecological assessment is a powerful, rapid and effective assessment strategy that enables complex spatio-temporal studies of community diversity and ecosystem function previously unavailable by traditional means.


Introduction
Modern human development has drastically increased the speed at which we alter our physical and societal environments, which can have quick and drastic effects on our ecosystems and their function.
We are often too slow to act on key changes in our environment, which makes the recovery and rehabilitation of healthy ecosystems more costly compared to programs that actively monitor ecosystems 1 . Active ecosystem monitoring relies greatly on the change in living biological communities (e.g. biodiversity) to assess changes in ecosystem function through community composition and ecosystem health, via community diversity 2,3 . However, despite the growing call to safeguard our natural ecosystems, we are currently experiencing a major decline in global biodiversity due to climate change and landuse alterations, which are outpacing current capabilities to actively monitor and respond to these threats 4 . It is therefore paramount that we develop more effective biodiversity assessment practices to increase our understanding of complex ecological systems and to promote ecosystem function and health 5 .
Accurate assessment of biodiversity relies on our understanding of the localized, spatial and temporal processes that shape changes in biodiversity in time and space [6][7][8] . Localized (i.e. site-speci c) biodiversity assessments largely predominate current monitoring practices, whereby local environmental conditions (e.g. landuse) are relied on to explain changes in biodiversity richness. As changes in biodiversity are also in uenced by spatial (e.g. dispersal) and temporal (e.g. phenology) factors, however, it is also essential to assess temporal and spatial community dynamics 6 . Furthermore, localized environmental ltering linked to biodiversity, which gives an indication of the general health of the ecosystem, can also be linked to ecosystem function if functional diversity of the localized biodiversity is available 2 . Biodiversity itself, can be quanti ed in many different ways, including through assessment of community richness or functional diversity. Richness is de ned as the number of unique taxonomic units per site/sample. Richness is often positively correlated with heterogeneous environments, which is often attributed to higher levels of functional diversity [9][10][11] . On the other hand, functional diversity directly quanti es the functionally disparate taxa within a community, and is becoming increasingly recognized as an important component of effective biomonitoring 12 . Additionally, differences in biodiversity (e.g. richness) between communities (e.g. sampling locations), commonly referred to as beta-diversity, determines if changes in biodiversity are in uenced by more local or spatial factors 13 . Partitioning the variance of beta-diversity into recently derived nested and turnover components provides even greater insight into the processes that are driving inter-community homogenization or dissimilarity 14 . It is therefore important to also assess changes in biodiversity at the local level and between communities to monitor changes in ecosystem health, and to assess functional diversity when able to determine the status of the local and regional ecosystem function. Unfortunately, current monitoring methodologies and practices are often limited or forced to simplify methods to cope with limited computational power or to reduce cost 15 . To implement increased spatial and temporal biodiversity assessment we have to develop and utilize improved biodiversity assessment methodologies, as current practices are limited in their e ciency to generate the data needed to rapidly assess ecosystems, particularly at increased spatial and temporal resolution.
The application of environmental DNA (eDNA) and metabarcoding has recently shown to offer increased sampling resolution for biodiversity assessment efforts [16][17][18] , though for some species, speci c traditional methods may outperform eDNA surveys 19 . Environmental DNA is extracted directly from an environmental sample (for example, water, soil or air) without prior isolation of the organisms themselves 20,21 . Sources of eDNA include sloughed skin cells, urine, feces, saliva or other bodily secretions, and could consist of both free molecules (extracellular DNA) and free cells [22][23][24][25] . Furthermore, eDNA collected from water samples has highly sensitive detection capability, is non-invasive to the sampled biota and not limited to physical environmental conditions (e.g. surface area or substrate types), thereby providing a wider sampling application with lower environmental impact, compared to traditional methods. It is important to note that eDNA capture can differ between environment types 26,27 , richness 27 and within sites 28 , due to seasonal variation, ecology, or random sample variation. Variability in e ciencies of traditional sampling methods between locations is well known, however, as the number of eDNA studies increase, so does our ability to account for random variability in data, particularly compared to early eDNA-based studies 18 . Combined with high throughput sequencing (HTS) applications, such as metabarcoding, eDNA based sampling is rapidly being integrated into standard ecological monitoring practices. Previous eDNA biomonitoring studies in natural environments have assessed population and communities across spatial and temporal scales for rivers [29][30][31] , lakes 28,32,33 and marine environments [34][35][36] . What is currently lacking is the link between functional and community diversity dynamics using eDNA across appropriate spatial and temporal scales. To develop and validate an eDNA-based approach Page 4/22 to biodiversity assessment it is therefore important to develop hypotheses from our current understanding of functional and community principles and dynamics.
Headwater riverine biodiversity is one of the longest standing realms of ecology and a key component of current freshwater biomonitoring and assessment 37 . Riverine catchments are a crucial component of regional biodiversity that harbor high levels of diversity due to their hierarchical structure, environmentally diverse habitats, and resulting unique headwater communities [37][38][39] . Freshwater macroinvertebrates are an invaluable source of ecosystem assessment information as their community assembly (localized) are strongly linked to ecological conditions 37,40 . Additionally, changes in macroinvertebrate communities in time and space also provides important information for any changes to ecosystems, resulting from effects of landuse alteration, environmental, or other anthropogenic effects 41 . Current biomonitoring practices in rivers utilize biological indices derived from freshwater macroinvertebrates 42 . Freshwater macroinvertebrates represent a wide range of species and functional groups, which respond dynamically to temporal and spatial environmental ltering, thus providing a clear depiction of localized ecosystem function 39,43,44 . Speci cally, functional feedings groups of freshwater macroinvertebrates provides direct assessment of nutrient cycling, productivity and decomposition 44 . Traditional taxonomic identi cation of freshwater macroinvertebrates, however, is largely limited to mature life stages that and can be di cult to identify or differentiate among similar species or genera 42 . The high level of taxonomic specialization required to identify specimens and the long processing times per traditional sample renders large scale ecosystem-wide traditional assessments expensive and time consuming 15,45 . In response there is a currently an ongoing rapid push to implement eDNA based riverine biodiversity assessment practices 21,26,30 , which forms the basis for this study.
Here, we assess seasonal patterns of biodiversity and functional diversity, using an experimental design which utilizes headwater sampling sites to associate local environmental conditions located within the same environmentally heterogeneous geographic region (e.g. catchment). We utilized a combined eDNA and traditional based biodiversity assessment approach to allow for direct comparison of historically supported ecological expectations from the traditional methods with molecular based eDNA methods. We investigated four main objectives and hypotheses. One, eDNA biodiversity will be greater and more diverse compared to traditional sampling, following previous experimental ndings. Two, riverine macroinvertebrate biodiversity, speci cally localized community richness, is expected to peak during spring and summer, when many stream macroinvertebrates are emerging as adults and reproducing, compared to fall and winter months when total biomass of many species has declined 28,46 . Three, utilizing the nested and turnover components of inter-community similarity (i.e. beta-diversity), we can expect high turnover within sites as community assembly changes over time, and high nestedness across environmental sites, attributed to environmental ltering of the localized sites. Alternatively, low nestedness could indicate a low effect of environmental ltering and a greater effect of stochastic or biotic factors in uencing the localized community assemblies. Four, environmental ltering effects that are linked to habitat modi cation, such as agriculture or urbanized areas, are expected to negatively impact macroinvertebrate diversity, particularly Chironomidae, and Ephemeroptera, Plecoptera and Trichoptera (EPT) taxa 47 , indicating variable site-speci c ecosystem health across the region. Functional diversity is also expected to be locally de ned with functional feeding groups being predominately shredders in undisturbed forested sites, and collectors in the sites dominated by ne particulate input, including urbanized and moorland habitats 48 .
We found eDNA biodiversity to be a better descriptor of the total macroinvertebrate diversity across all sampling sites, with trends between the two methods showing general similarities. Community richness was greater in the spring and summer and lowest during the winter, as expected with the greatest change in community composition across seasons linked to changes in Chironomidae genera richness. Landuse, while showing distinct environmental differentiation, was not associated with local community richness.
Additionally, spatio-temporal dynamics among communities were found to be predominately turnover driven, indicating strong seasonal or region-wide effects, whereas nestedness effects were mostly limited, suggesting weak localized environmental sorting of the communities. Lastly, functional diversity showed clear region-wide generalization of feeding functionality, suggesting biodiversity is driven by regionalbased bottom-up dynamics, which suggest biodiversity management should focus on regional over localized spatial extents.

Bioinformatics
After stringent ltering and quality control, 12,592,362 reads were obtained with an average of 74,954 (± 31,050) reads per sample. Negative controls, while showing no bands on agarose gels post library preparation, generated 676 reads across all blanks (N=12). Of the negative control reads, 411 were unknown bacteria, 3 reads were associated with three genera of Rhodophyta (red algae), 2 reads were linked to unknown fungi, and 260 reads were linked to a single Dipteran ASV across four blanks. For downstream analyses, the Dipteran ASV was removed from subsequent analyses, whereas all other potential contaminants were not included as they were non-targeted. In total, 20,437 ASVs were identi ed. Average reads per site, after rarefaction was 75,871 (± 37,670) with four sites having less than 10 000 reads from four different sites across three different landuse types and two seasons. The average number of taxonomic assignments was 13,260 (± 9,567). The average number of kicknet sampling specimens per site was 1,529 (± 1 555). Mean site singletons were 9.86 (± 6.67) for eDNA and 9.02 (± 3.29) for kicknet sampling. Singletons were included in subsequent analyses given the robust use of sample replication used in the study design, whereby if a sequence was not observed in at least 2 of the 3 replicate samples the sequence was not included in the downstream analysis.

Community dynamics
All environmental variables and their associated summary statistics are presented in Table 1 & Figure 2.
Overall, we observed 226 unique genera using the eDNA based approach and 83 genera using the traditional kick-netting approach ( Table 2). On average, eDNA genera accounted for 78.2% of the unique observed diversity in a given site, with traditional methods accounting for 5.9% with an overlap between the two methods of 15.9% ( Figure 3). Key differences between the methods were the higher number of genera observed using eDNA vs traditional methods for Chironomidae (75 vs 10), Oligochaeta (23 vs 2), Trichoptera (33 vs 24), Rotifera (8 vs 0), Coleoptera (20 vs 14) and Copepoda (5 vs 0) ( Table 2). The full breakdown of genera per landuse type per sampling method can be found in Supplementary Table S1.
Genera richness derived from eDNA was signi cantly greater than traditional derived eDNA across season (p<0.01) and landuse type (p<0.001). There was a signi cant landuse x method interaction effect (p<0.01) and season x method effect indicating non-covarying biodiversity dynamics between the two methods. Results for EPT also indicated non-covarying biodiversity dynamics with signi cant landuse x method (p<0.001), and method x season (p<0.001), whereas both methods showed signi cant differences across seasons and landuses. Chironomidae diversity also showed signi cant landuse x method (p<0.001) and method x season (p<0.001) interactions as well as signi cant effects of landuse and seasons on richness dynamics. Functional diversity showed signi cant landuse x method (p<0.001) and method x season (p<0.001) interactions, as well as seasonal (p<0.001) ( Table 3 & Fig. 4).
Change in community and functional diversity over time and space Turnover was signi cantly greater than nestedness across landuse (p<0.001) and season (p<0.001) for genera derived from or traditional methods ( Figure 5). The variation between turnover and nestedness did not differ across season (p=0.66) or landuse (p=0.082), but did differ signi cantly between methods (p<0.001), with eDNA showing greater greater sensitivity to detect nestedness versus traditional methods ( Figure 5). Positive and negative relationships with regards to landuse relate to the environmental gradient depicted in Fig S1 and described above. Changes in biodiversity over time, here turnover dominated, were directly related to the change in functionality over time (Fig. 6). For eDNA samples the increase in overall functionality was largely driven by the increase in number of genera from Diptera (all landuses), Coleoptera (Acid grasslands, Moorlands, Urban and Agriculture), Ephemeroptera (Forest), and Trichoptera (Forest and Moorlands). Transitioning into summer, acid grasslands and moorlands showed a loss of scraper and collector functionality stemming from losses in Plecoptera and Trichoptera genera. Losses of Plecoptera genera in forest sites resulted in a loss of scraper functionality, whereas gains in Plecoptera in the urban environments showed gains in collector functionality. Agricultural sites showed gains in scraper and gatherer function, due to some increases in Diptera and Trichoptera. Transitioning into fall indicated loss in functionality across landuse, driven by losses in Diptera (all landuses), Coleoptera (all landuses), Ephemeroptera (agricultural, moorland), and Trichoptera (urban, forest, agricultural). Winter was largely static, with the exception of gains in functionality for agricultural sites, which was driven by increased occurrence of Plecoptera and Trichoptera genera. Traditional samples indicated slight increases in functionality for agricultural and moorland landuses, primarily from increases in Plecoptera (agricultural, urban) and Coleoptera (moorland) genera. Summer increases in functionality for acid grasslands and moorlands were linked to increases in Trichoptera and Plecoptera in acid grasslands, and Plecoptera and Coleoptera in moorlands. Declines in fall gatherer function, in agricultural and forest landuses, were driven by losses in Trichoptera (agriculture sites) and Plecoptera (forest sites). Loss of functionality in winter for acid grassland, agricultural and moorlands was driven by losses in Trichoptera (acid grasslands), Plecoptera (acid grasslands, moorlands), Ephemeroptera (acid grasslands, agricultural, moorlands), and Coleoptera (acid grasslands).

Discussion
We show that eDNA based assessment offered a ner resolution of the spatial and temporal dynamics, making it a more e cient means to assess biodiversity dynamics. Additionally, biodiversity patterns derived from eDNA and traditional methods showed similar temporal trends in richness and functional diversity. Localized biodiversity richness showed no environmental ltering across the landuse types, whereas partitioning of beta-diversity showed clear differences in spatio-temporal dynamics. Speci cally, regional environmental conditions were the main driver of biodiversity change and landuse effects were less pronounced. Importantly, we show here that eDNA based biodiversity assessments show meaningful spatial and temporal relationships, including increased ability to detect important indicator taxa, particularly Chironomidae and EPT taxa, which are di cult to directly sample for many of the locations or at different time points in the year. We also show the rst example of eDNA derived functional temporal spatial dynamics, which provides clear management related information on how the ecosystem can be effectively managed.
As expected, biodiversity was greater during spring and summer months for both sampling methods. Genera richness dynamics differed between eDNA and traditional methods among landuse types, but not among seasons where eDNA derived diversity was consistently greater than traditional diversity. Importantly, we found greater eDNA derived diversity compared to traditional methods for all sites. While several studies have shown increased observable biodiversity using eDNA over traditional methods 17,32 , they have predominately focused on sh, whereas macroinvertebrate focused studies vary slightly with most showing greater 30,49 , but also some with lower eDNA diversity 50 , compared to traditional sampling. The disparity in macroinvertebrate diversity may stem from the increased di culty in designing a suitable primer to capture the full range of diversity associated with macroinvertebrates, though see Leese et al. 2020 51 . Additionally, across landuse types, acid grasslands and moorland sites had greater eDNA diversity compared to agriculture, urban or forest sites. The results from eDNA monitoring were more in agreement with common expectation that unmodi ed landuse should hold greater biological diversity, particularly with regards to the higher diversity found in moorland and acid grassland sites, which were the least modi ed areas in the catchment. Traditional sampling, however, suggested greater biodiversity in agricultural and urban sites compared to moorlands or acid grasslands. The lower biodiversity observed with traditional methods metrics in the less disturbed sites is likely due, in part, to substrate types. The moorland and acid grassland sites are dominated by large boulders or loose sediment, which are not ideal substrates when performing traditional kicknet sampling methods, that perform better with gravel substrate (signi cant positive richness with increased gravel coverage p = 0.007), and may lead to under sampling of the local taxa 45 . Previous studies have suggested eDNA transport from upstream communities can increase sampled biodiversity 52,53 . There is, however, no clear indication that the landscape or catchment surrounding these sites is any more diverse, as the moorland sites in particular are at higher elevation and more isolated compared to the other landuse sites. With upstream transported limited across the study system, the increased biodiversity is more likely an effect of eDNA versus traditional methodology than proposed eDNA ecology. As a recommendation, eDNA sampling is likely to be more bene cial overall compared to traditional sampling for detecting higher biological diversity, which has direct implications for detecting traditionally harder to identify biomonitoring groups such as Chironomidae or Diptera. Additionally, the use of eDNA is highly bene cial for sampling in traditionally di cult to sample locations, including non-traditional substrate types which could introduce bias when traditional, or bulk, sampling methodologies are used.
Environmental DNA biodiversity was greater for general biodiversity and all other subsets of biodiversity, including EPT, chironomids and functional diversity (Figures 3 & 4). The greatest increase in eDNA biodiversity resolution was in traditionally hard to identify groups, which would otherwise not be observable using taxonomic derived methods. Most revealingly, eDNA data more accurately depict Chironomidae life cycle patterns compared to traditional sampling. Speci cally, seasonal variation in eDNA derived richness follows the expected larval emergence patterns, which increase over spring and summer, and steadily decline over fall and winter 28 , in contrast to traditional sampling which was unable to detect this seasonal shift, possibly due to the di culty in traditionally observing Chironomidae. Likewise, EPT emergence was detectable via both eDNA and traditional sampling, as Ephemeroptera generally emerge during the spring whilst Trichoptera emerge at various times throughout the year.
Spatio-temporal dynamics showed biodiversity was driven by regional turnover dynamics and less affected by localized environmental specialization or nestedness (Figures 5 & 6). Nestedness across sites was not signi cant with either method, suggesting that differences in biodiversity between sites was not due to localized environmental sorting, as per our initial expectations. The primary driver of observed heterogeneity in between site biodiversity was largely due to the seasonal turnover, likely driven by the high disturbance events historically occurring in the region of the study 54 , or due to very stronger effects of biotic interactions 8 . For the eDNA based methodology, there was a greater turnover observed in urban and forest sites, which was attributed to higher pH and lower moss and boulder coverage compared to other sites. Turnover in EPT and Chironomidae showed similar trends compared to overall turnover, whereby the differences in biodiversity between sites was signi cantly attributed to seasonal replacement of genera along the environmental gradient, which was predominately linked to pH and substrate type ( Fig S1). Conversely, turnover was greater in moorland sites, linked to increased boulder and moss coverage. The disparity in the observed relationship between methods is likely driven largely by the methods themselves and the underlying richness values for each method, as mentioned above. Overall, both traditional and eDNA based turnover suggest seasonal turnover dominated differences in biodiversity, which could be attributed to adaptation (historical) of communities to regional conditions 55 . The biodiversity across the system is more likely a product of frequent founder and colonizing effects resulting from frequent disturbance patterns, and inability of the sites to establish long-term interacting communities 56,57 .
The functional diversity was largely dominated by collector feeders, indicating that the regional biodiversity assembly is driven by ne particulate organic matter (FPOM), which is also referred to as seston 58,59 . Two main factors contribute to the high FPOM driver across the study area. First, the widespread fecal input from animal agriculture that covers the entirety of the catchment. Additionally, for upland sites where agriculture is less prevalent, moorlands produce a high amount of FPOM 60 . The combined FPOM inputs from moorland and agricultural/urbanized environments creates a regional-wide FPOM system, thereby homogenizing the functional habitat of the region as a whole. Whereas the moorland effect is likely limited to certain headwater sites in this study, the quality difference in FPOM generated from agricultural vs moorland sources likely plays a role in the diversity differences seen between sites where collectors dominate on the whole 43,58 . Seasonal shifts in functional traits were observed with eDNA, but not with traditional methods. This is re ected in the ability of eDNA to detect more functional groups compared to traditional methods, particularly with regards to Diptera, including Chironomidae and Simuliidae (i.e. black ies), which are important collector groups that are strongly affected by changes in temperature 60 . Likewise, eDNA functional assessment indicates a change in grazer functionality with season, closely following expectations of periphyton availability, which is a key driver of grazer activity 61 . The homogeneity in functional diversity between the landuse sites for both eDNA and traditional based methods further indicated that environmental ltering was not the primary driver of biodiversity difference between sites. A key take home message from assessing functional diversity in the study, over simply relying on variation in richness, is that these ndings point to a clear management strategy to increase diversity across the system. Speci cally, regional biodiversity would bene t by increasing the habitat for collector functional groups through improved management of local broadleaf forest, and agricultural practices to increase coarse feeding material at key head water sites, which would increase the overall ecosystem stability of the region by increasing environmental heterogeneity.
A wider implication arising from this study is the suggestion that eDNA can disclose a much greater resolution of diversity compared to traditional approaches and can enable multiple levels of analyses from a single data set. The bene t of comparing traditional ndings with eDNA based analyses, which can be analyzed using standardized approaches, is immensely valuable and will help avoid unintended biases introduced from cross study traditional protocols. Overall, our ndings show that eDNA is a more effective survey method to sample macroinvertebrates and provides clearer indications of the seasonal and environmental effects on multiple levels of diversity compared to traditional methods. Additionally, we provide a key assessment of regional biodiversity dynamics, which are currently underrepresented in the literature. Speci cally, we show that the increased resolution of eDNA based biodiversity assessment effectively separates spatio-temporal and localized biodiversity dynamics, which here shows the importance of regional over localized management strategies. Determining such regional drivers allows for effective biological management, whereby ood control, versus altering current landuse practices, is more likely to have a greater impact on biodiversity in disturbance driven sites. Finally, by utilizing functional diversity assessment we show a clear reason for why community composition has arisen. Empowered by eDNA metabarcoding and appropriate ecological synthesis, we therefore provide a more valuable means to describing biodiversity than simply counting unique individual units with no link to what the numbers mean at ecological and ecosystem scale.

Study Area
The Conwy Catchment is a 678 km 2 river drainage in north Wales that encompasses a wide range of habitats including forest, moorland, agriculture, light urbanization, and acid grasslands (Fig. 1). The area experiences rapid climatic shifts, particularly during winter months due to its mountainous terrain and porous rock foundations, which facilitates ash ooding, making it susceptible to periodic disturbance 54 . The Conwy Catchment area exhibits four distinct seasons that correspond to the expected life cycles of EPT and Chironomidae larval emergences throughout the year.

Sampling
Fourteen headwater sites across 5 landuse types (acid grassland, agriculture, forest, urban, moorland) were sampled once per season (spring, summer, fall, winter), during 2017 (5 landuse types x14 sites x4 seasons). Headwater sites were selected to ensure local landuse types were not in uenced by other landuse type via downstream transport 62,63 . In total 168 eDNA samples and 56 traditional kick-net samples were taken over the course of the study. Sampling for each season occurred over two consecutive days. For all sampling events, streams were sampled for both eDNA and macroinvertebrate community composition during the same day. Water samples for eDNA analysis (1 L) were collected in triplicates from each stream with plastic bottles that had been cleaned prior using a 10% bleach solution (soaked for 1 hour). These were then ltered through 0.22 μm SterivexTM lter units (EMD Millipore Corporation, Billerica, USA) using a Geopump TM Series II peristaltic pump (Geotech, Denver, USA). As lters would occasionally experience reduced ltration e ciency for different sites, or seasons, due to stream sediment loading, we would continue to run the pump for each sample until at least 500ml was ltered, to avoid potential downstream variation in sampling 64 , which we previously showed was not an issue in this experimental setup 53 . The lters were immediately preserved in 500 μl ATL lysis buffer (Qiagen, Venlo, The Netherlands) and stored in coolers during same-day transit to the laboratory, then stored at 0 °C, for further processing. Macroinvertebrate communities were sampled using a standardized 3-minute kick sampling protocol, with a 500 μm mesh gauge kick/hand net. Kicknet sampling occurred after eDNA sampling to ensure disturbance of the site from kicknetting would not in uence the eDNA signal. Both bank margins and ri e habitats were sampled during this timed sampling period. Macroinvertebrates were preserved in absolute ethanol (99.8%; VWR International, Lutterworth, UK) on collection. Upon return to the laboratory the macroinvertebrates were cleared of other collected material and identi ed to the lowest practical taxonomic per protocol 65 . Environmental (abiotic) data including pH, conductivity, depth, moss coverage, algal coverage, plant coverage boulder coverage, gravel coverage and sand coverage, were collected for each site following UK Environment Agency site assessment protocols 65 . Three eld blanks consisting of deionized water were prepared in Nalgene sampling bottles and were kept with the other sampling gear throughout each sampling event. For each seasonal sampling event we incorporated three eld blanks into the sampling process, processing one blank every 4 th sampling site, resulting in a total of 12 blanks for the study. Field blanks were 1L in volume and treated the same as standard samples.

Extraction and Sequencing
We followed unidirectional lab practices from eld, to extraction to library preparation by using designated extraction (PCR free) and library preparation rooms. DNA was extracted from the lters using a modi ed QIAGEN DNA blood and tissue extraction protocol 66 . In short, 70 μl proteinase K was added directly to the lters and incubated at 58 °C overnight in a rotating hybridization chamber. Then, the lysate was extracted and the full volume was ltered through a spin column tube, after which point the standard extraction protocol was continued. Extracts ( nal volume 50 μl) were then cleaned for impurities using QIAGEN Power Clean kit and frozen at -20°C for subsequent analyses. Sequencing libraries were created using a two-step protocol (see Bista et al. 2017), using matching dual end index tags (IDT) and the following COI gene region primers for the rst round of PCR (PCR1): m1COIintF (5′-GGWACWGGWTGAACWGTWTAYCCYCC-3′) and jgHCO2198 (5′-TAIACYTCIGGRTGICCRAARAAYCA-3′) 67 . Libraries were created at Bangor with the assistance of a Gilson pipette max liquid handler before being shipped to University of Birmingham's Genomic sequencing facility for quality control and sequencing.
Round 1 ampli cation (PCR1) with the COI primers was performed in triplicates, which were then cleaned for primer dimers using magnetic beads (Beckman coulter), pooled and index labeled during Round 2 PCR step (PCR2), which were cleaned again using magnetic beads. Unique dual paired end-indices were designed and purchased from Integrated DNA Technologies, to complement the Illumina P5/P7 sequence adapters. PCR1 utilized Thermo Scienti c's Ampli-gold mastermix due to the high number of inosine in the COI primer pair, and for PCR2 we utilized New England Biolab's Q5 mastermix. All PCR1 and PCR2 reactions were run in 25 μL volumes. PCR1 amplicons were generated using a reaction mix of 12.5 μL mastermix, 2 μL DNA template, 1 μL of each primer and 8 μL nuclease free water and ampli ed using an initial 95 °C for 5 min then 25 cycles of 95 °C for 30 s, 54 °C for 30 s and 72 °C for 60 s followed by a 72°C nal annealing for 10 min. PCR2 amplicons were generated using a reaction mix of 12.5 μL mastermix, 2 μL DNA template, 1 μL of each primer and 8 μL nuclease free water and ampli ed using an initial 98 °C for 30 s then 15 cycles of 98 °C for 10 s, 55 °C for 30 s and 72 °C for 30 s followed by a 72 °C nal annealing for 10 min. The, PCR2 amplicons were puri ed using High Prep PCR magnetic beads (Auto Q Biosciences) and quanti ed using a 200 pro plate reader (TECAN) with the Qubit dsDNA HS kit (Invitrogen). The nal amplicons were pooled in equimolar quantities (at a nal concentration of 12 pmol) using a Biomek FXp liquid handling robot (Beckman Coulter). Pool molarity was con rmed using a HS D1000 Tapestation ScreenTape (Agilent). Sequencing was performed on an Illumina HiSeq platform 250bp Paired-End, with an intended coverage of 100,000 reads per sample.

Bioinformatics
Bioinformatic processing up to taxonomic assignment was performed by University of Birmingham. In short, per base quality trimming was performed on demultiplexed reads using SolexaQA++ v. 3 69 allowing up to 3 mismatches per primer sequence. Only sequences with both forward and reverse primers were retained for further analyses. Amplicon sequence variants (ASVs) were obtained via Usearch at 97% similarity threshold, and denoising with the -unoise3 algorithm. Chimeras were removed as part of the -unoise3 algorithm 70 .
Taxonomy to the genus level was assigned to representative ASV with BLAST against ASVs using the non-redundant nucleotide database of NCBI, using the default settings 71 .

Metacommunity analyses
All statistical analyses were performed using R version 3.6.1 72 . Sequence reads were rari ed for each set of replicates to the lowest replicate level. Mean number of reads for each ASV were calculated across the sample replicates before being matched to their taxonomic identi er. ASVs that were not identi able to the genus level, or to a functional group (below), were not included in subsequent analyses. Genera richness was calculated as the number of unique genera per site. We further divided richness into unique EPT genera for EPT richness, and Chironomidae richness as the number of unique Chironomidae genera. Functional richness was calculated as the partition of unique functional groups per sample, following the partition of functional groups in Moog (2017). In short, Moog (2017) provides a catalogue of 3 296 metazoa species that form the basis of ecological status assessment for many European environmental agencies. The functional scores are assigned to each taxa based on a ten point partitioning to re ect the variation in functionality within taxa, meaning the functional scores is a score of the function and not simply a reassignment of the taxonomic identi cation. These groups re ect the functional feeding groups, divided into 8 categories, including; grazer/scrapers, xylophagous, shredders, gatherers/collectors, active lter feeders, passive lter feeders, predators and other. We re ned these groups to shredders, grazers, and collectors whereby collectors were the summation of gather/collectors and lter feeding groups to simplify the functional groups to those used across wider studies. To assess community dynamics between sites we calculated the nestedness and turnover components of betadiversity following Baselga (2010) 14 , whereby beta-diversity was calculated as βsor (1), Turnover as βsim (2) and βnes as βsor minus βsim (3) 14 . All mathematical formula follow the nomenclature of Baselga (2010), with a being the number of genera common to both sites, b is the number of genera occurring in the rst site but not the second and c is the number of genera occurring in the second site but not the rst.
We used generalized least squares (gls), as implemented using the gls function in the nlme package 73 , to assess the statistical relationships between (1) community, (2) EPT, (3) Chironomidae, and (4) functional richness (each as a separate response variable and independent statistical test) and all two-way interactions of the explanatory variables, including sampling method (eDNA or traditional), season (spring, summer, fall and winter) and landuse gradient (see below for description). Generalized least squares is an extension of linear regression that allows for variance structuring of variables to account for suspected correlation between residuals. Here, we speci cally used the gls framework to account for potential spatial autocorrelation between sites by included a variance structure using the distance matrix for the sampling sites 73 . We further used gls to assess the statistical relationships between nestedness and turnover between communities against a set of explanatory variables, including the landuse gradient, sampling method, and season, including all possible two-way interactions. Backward model selection was performed to nd the most parsimonious model using Akaike information criterion (AIC) to determine the best model t 74

Data availability
Data including amplicon sequence data, ASV matrices and environmental data associated with the manuscript will be deposited on gshare.
Code availability