Effects of fungicides on aquatic fungi and bacteria: a comparison of morphological and molecular approaches from a microcosm experiment

Fungicides are frequently used in agriculture and can enter freshwater ecosystems through multiple pathways. The negative impacts of fungicides on microorganisms, fungi in particular, and their functions such as leaf decomposition have been repeatedly shown. In our previous microcosm experiment with three consecutive cycles of fungicide exposure and colonisation of leaf substrate, we found clear functional changes, but no differences in fungal community structure could be detected using morphological identification by analysing the spores of aquatic hyphomycetes. In this study, we examined the effects on fungal and bacterial community composition in detail using ITS and 16S metabarcoding and comparing the results to morphologically assessed community composition. While we found fewer species with metabarcoding than with morphological identification, metabarcoding also enabled the identification of several fungal species that were otherwise unidentifiable morphologically. Moreover, by distinguishing individual amplicon sequence variants (ASVs) metabarcoding provided greater taxonomic resolution. In line with the morphological results, metabarcoding neither revealed effects of fungicides on the aquatic hyphomycetes nor on the total fungal or bacterial community composition. However, several ASVs responded significantly to fungicides, demonstrating variable tolerances within species. Overall, the absence of detectable effects of fungicides on the community structure despite clear functional effects, suggests a complex relationship between community structure and the ecosystem function of leaf decomposition.


Background
Freshwater ecosystems such as streams host a disproportionately large part of global biodiversity [1].They are strongly influenced by anthropogenic stressors, including habitat degradation and pollution, which are typically related to land use within their upstream catchments [2][3][4].One major group of pollutants are pesticides originating from diffuse sources such as spray drift, drainage and surface runoff from agricultural land use [5][6][7] as well as from point sources, such as wastewater treatment plant discharge (from urban and agricultural land use; Müller et al. [8]).Among pesticides, fungicides are broadly applied in agriculture against a wide range of fungal pathogens and they widely occur in aquatic systems [9].Since fungicides target basic biological processes in organisms, such as energy production, they pose a risk to a broad range of aquatic life, especially fungi [9].Moreover, fungicide exposure is predicted to increase with global warming and the related spread of fungal diseases as well as the development of fungicide resistance [10,11].Yet, major research gaps remain regarding the effects of fungicides on microbial communities and their associated function, especially as interactions between microbial taxa can affect decomposition [12].
Decomposition of allochthonous organic matter, such as leaves, is a crucial ecosystem process of energy provisioning in headwater streams given that shading often limits primary production [13,14].Leaf decomposition is controlled by abiotic (e.g., temperature, flow velocity, sediment, dissolved chemicals and other pollutants) and biotic (e.g., leaf characteristics, microbial and invertebrate communities) factors [15][16][17].Fungi and bacteria are the main microbial groups contributing to decomposition.Typically, fungi contribute more to leaf decomposition than bacteria and dominate the microbial biomass associated with leaves [18][19][20].The saprotrophic fungal communities in streams are often dominated by aquatic hyphomycetes [21], which are particularly important for decomposition [22][23][24].Hyphomycetes are a diverse, polyphyletic group with a cosmopolitan distribution [25] and are characterised by large, complex, asexual spores (conidia).In contrast to studies based on classical methods, such as morphological identification of such spores, molecular studies consistently demonstrated the ubiquitous presence of other fungal taxa [26,27], whose contribution to decomposition are still largely unknown.Therefore, a better understanding of fungal diversity (the 'mycobiome') and the relative contribution of different fungal taxa to leaf decomposition are required to predict the consequences of anthropogenic stressors for microbial leaf decomposition.Furthermore, as bacteria interact with fungi, impacts on fungi, e.g., through fungicide exposure, can indirectly affect bacteria [12,28].
The classical morphological identification of spores of aquatic hyphomycetes is limited by some taxa not growing or releasing spores, while the identification of other taxa remains ambiguous due to a lack of diagnostic characteristics or phenotypic plasticity [29,30].Therefore, these taxa, while potentially important in the fungal community and associated ecosystem functions, may be missed via morphological identification [26].A modern alternative is a taxonomic identification based on genetic markers, which is independent of the presence of reproductive structures, such as spores.Such identification methods rely on polymorphic gene regions to distinguish between different phylotypes (e.g., [31]) and ultimately species (DNA barcoding [32]), where a taxonomic assignment is performed through comparisons to reference databases [33].Genetic markers may provide better resolution and allow the identification of taxa that are missed in morphological identification due to low or no sporulation or lack of identification characteristics [34,35].The amount of unidentified fungal taxa retrieved via DNA barcoding of whole communities in a sample (i.e., DNA metabarcoding [36]) can be considerably higher than sporulation-based assessments [37].
For the taxonomic assignment in metabarcoding studies, comprehensive and well-curated reference databases are essential.For fungi, different genes of the rRNA operon were used in the past, leading to several reference databases with only partly overlapping species repositories [38].Thus, DNA metabarcoding studies using only one short fragment (< 500 bp) from among the possible barcoding gene fragments can underestimate the known species diversity when assigning Molecular Operational Taxonomic Units (MOTUs) [39].Long-read sequencing is one of the approaches devised to mitigate this problem by linking the disparate information of different databases created from different fragments and thereby improving the taxonomic assignment [39].Besides the assignment to species or strains with defined thresholds (= MOTUs, threshold often 97% sequence similarity), metabarcoding allows exploring intraspecific diversity using Amplicon Sequences Variants (ASVs) based on metabarcoding data [40].
In a recent microcosm study, we found a significant change in microbial leaf decomposition under fungicide exposure, which was not reflected in changes in community composition, where the community was assessed via morphological identification of aquatic hyphomycetes [41].In addition, sporulation data of this study was incomplete due to device malfunctioning.This inhibited understanding of the mechanisms that may govern the association between functional and structural change, i.e., changes in decomposition and community composition.In the present study, we, therefore, used DNA metabarcoding with ITS and 16S as genetic markers to further examine the response of leaf-associated fungal and bacterial decomposer communities to multiple cycles of fungicide exposure and resource colonisation, investigating the structural and functional responses and unravelling potential associations.We hypothesised that: 1. using metabarcoding, most morphologically identified fungal taxa will also be identified, but a much greater diversity as well as overlooked aquatic hyphomycete taxa and other fungal taxa will be identified; 2. by accounting for overlooked fungal taxa in the fungal community, significant structural and functional responses of the fungal community to repeated fungicide exposure will be detected; 3. similar to hypothesis 2, significant responses of the bacterial community and its functional composition to repeated fungicide exposure will be detected using metabarcoding.

Experimental design
To study the response of leaf-associated fungal and bacterial decomposer communities to repeated fungicide exposure, we conducted a microcosm experiment (see [41] for a detailed description, Additional file 1: Fig. S1).We used 14 microcosms, each containing 4 L of pre-filtered water from an unpolluted stream (Sulzbach, Rhineland-Palatinate, Germany).To study decomposition, bags containing 4 g (± 0.01 g) of black alder (Alnus glutinosa (L.) Gaertn.)leaves in nylon mesh bags (mesh size: 2 mm) were placed in a forested stream without agriculture in the upstream catchment for 7 days to enable microbial colonisation by local communities (bags of cycle 1) and then transferred to the microcosms.For cycles 2 and 3, the leaf bags (mesh size 10 mm) were directly inserted into the microcosms for 7 days to allow for colonisation by the microbial communities of the previous cycle.Following colonisation, a short-term 48-h fungicide peak consisting of a fungicide mixture of four fungicides (metalaxyl, prothioconazole, pyrimethanil, and prochloraz) detected in streams in the same geographic region was applied at a field-realistic level to half of the microcosms [9] (nominal fungicide peak concentrations between 0.2 and 40 µg/L), resulting in seven experimental replicates per treatment (for details see [41]).The peak was followed by 12 days of exposure to baseline levels of fungicides at tenfold lower concentrations [42].This resulted in a cycle length of 3 weeks, while the whole experiment lasted for 7 weeks (Additional file 1: Fig. S1).Five of the seven replicates per treatment (i.e., control and fungicide exposure) were sampled for morphological and DNA metabarcoding analyses after the end of each cycle, resulting in a total of 30 samples.Thirteen leaf disks with 1 cm diameter were randomly cut from the leaves avoiding the midrib within each bag sample at the end of each of the three cycles.The remaining material was dried to a constant weight at 60 °C, weighed to the nearest 0.01 g and the decomposed leaf mass per degree day was calculated for each replicate.For further details, see Schreiner et al. [41].

Determination of bacterial density
Bacterial cell density was determined using epifluorescence microscopy following the approach of Buesing [43].Briefly, bacteria were detached from five formalin-preserved leaf disks by ultrasonication, filtered on an aluminium oxide membrane filter (0.2 µm pore size, Anodisc, Whatman), stained with SYBRGreen II and the cell numbers were automatically counted using an imageanalysis program (AxioVision 4.8, Carl Zeiss).The resulting counts were normalised to the processed leaf area.

Morphological identification of aquatic hyphomycetes
Aquatic hyphomycetes were morphologically identified from sporulated conidia as part of a previous study [41].To induce the sporulation of aquatic hyphomycetes, the five leaf disks from each sample were submerged in water of the respective microcosm and orbitally shaken at 120 rpm in darkness at 16 °C (the same temperature as used during the experiment).After 72 h, conidia were prevented from agglomerating using 0.5% Tween80 and fixed in formaldehyde (final concentration 2%).The samples were then stored at 4 °C.Immediately before identification, an aliquot of the conidia suspension was homogeneously filtered through a gridded membrane filter (0.45 μm), and retained conidia were stained with lactophenol cotton blue.At least 300 conidia were identified for each replicate (100-400× magnification) primarily using the key by Gulis et al. [21].The number of counted conidia was normalised to the total filter surface, used sample volume and dry weight of the respective set of leaf disks.A malfunction of the orbital shaker resulted in minimal sporulation for the first cycle, and consequently, the related data are missing.

DNA extraction
For each sample, three leaf disks were stored at − 80 °C until DNA was extracted using the Qiagen DNeasy PowerPlant Pro Kit following the manufacturer's instructions (samples were homogenised by vortexing).Three microliters of extracted genomic DNA were checked on a 1% agarose gel, concentration quantified using a Qubit BR dsDNA Assay (Invitrogen) and then diluted to 10 ng/μL of DNA using molecular-grade water.DNA was also extracted from the 21 fungal isolates (Additional file 1: Table S1) using the Qiagen DNeasy PowerPlant Pro Kit (Additional file 1: Fig. S1).

Fungal ITS metabarcoding
For hyphomycetes, Letourneau et al. [35] reported the best interspecific separation using ITS compared to other genes; therefore, the ITS2 fITS7 [44] (5′ GTG ART CAT CGA ATC TTT G 3′) and ITS4 [45] (5′ TCC TCC GCT TAT TGA TAT GC 3′) primer pair was used to target a 122-245 bp fragment of the ITS2 region.The fITS7 primer also tends to be more fungi-specific, amplifying less plant DNA (see the review article by Lear et al. [46]).
To prepare the amplicon libraries, two rounds of PCR were performed: to amplify the target region and to attach sample-specific indexes as well as Illumina sequencing adaptors.For the first PCR, the forward primer fITS7 (ACG ACG CTC TTC CGA TCT GTG ART CAT CGA ATC TTT G) and the reverse primer ITS4 (CGT GTG CTC TTC CGA TCT TCC TCC GCT TAT TGA TAT GC) with overhang adapters/linkers attached were used in the following PCR reaction: 1 μL of DNA, 0.5 μL of each 10 μM primer, 0.5 μL 10 mM deoxyribonucleotide triphosphates (dNTPs), 2.5 μL 10 × 5 Prime HotMaster Taq Buffer with Mg 2+ and 0.25 μL 5 Prime HotMaster Taq DNA polymerase (5 U/µL, Quantabio) in a total volume of 25 μL.The reaction conditions were as follows: an initial denaturation step of 94 °C for 5 min, then 25 cycles of 94 °C (30 s), 52 °C (30 s) and 72 °C (45 s), and a final extension of 72 °C for 10 min.Two technical PCR replicates were performed.The success of the first PCR was determined on a 1% tris-borate-EDTA (TBE) agarose gel.The products of the first stage amplification (pooled technical replicates) were then purified using 1.8× Agencourt AMPure XP magnetic beads (Beckman Coulter) to remove free primers and primer dimers.The second PCR was performed with primers designed based on the Illumina TruSeq and Nextera PCR primers.This reaction was run similar to the first stage except ten cycles were used instead of 25.The success of the second PCR was again determined by running a 1% TBE agarose gel.The products of the second PCR were then cleaned for a second time using AMPure XP magnetic beads.Concentrations of the samples were then measured using a Qubit fluorometer (Thermo Fisher) and the libraries were pooled equimolarly in two batches (half the samples in each).Final libraries were quality checked with a Fragment Analyser (Agilent) before sequencing on two 2 × 250 bp MiSeq V2 runs at Eurofins GmbH (Constance, Germany).

Bacterial 16S metabarcoding
For bacteria, the V4 region of the 16S SSU rRNA bacterial gene was amplified using the 515F forward primer [47] and the 806R reverse primer [48].The 16S PCR amplification was conducted using a similar two-step protocol as described for ITS above, except for the conditions of the first PCR: initial degradation at 94 °C for 3 min; 25 cycles of 94 °C (45 s), 50 °C (60 s) and 72 °C (90 s); final extension at 72 °C for 10 min.The libraries were pooled equimolarly in two batches and sequenced on two 2 × 250 bp MiSeq V2 runs at Eurofins GmbH (Constance, Germany).

Long-read PCR and library preparation of fungal isolates
Strains of 21 fungal species (including 17 species of aquatic hyphomycetes) (Additional file 1: Table S1) known to occur in local streams, were used to generate reference sequences.Fungal strains were obtained from various streams in Germany as described in Baudy et al. [49], via isolation and cultivation of single spores (following the method described in Descals [50]).The strains were axenically cultured on malt extract agar (20 g/L malt extract, 20 g/L agar) at 16 °C in darkness and subcultures were deposited at the Leibniz Institute DSMZ (German Collection of Microorganisms and Cell Cultures GmbH, Braunschweig, Germany).

Sequence data processing
The 16S and ITS amplicon data were analysed using the DADA2 pipeline version 1.14.1 [40] and low-quality and chimeric sequences were removed following the steps described by Pauvet et al. [52].For the 16S data, forward and reverse reads were merged, but due to the low quality of the reverse ITS reads, the majority of the forward and reverse reads could not be merged and, therefore, only forward reads (R1), truncated to 150 bp, were taken for analysis.Taxonomy was assigned using the UNITE [33] and SILVA version 132 [53,54] databases.Additional filtering was carried out using the R package phyloseq version 1.30.0, to remove 16S rRNA gene sequences annotated as chloroplast (order level) or mitochondria (family level).Relative read abundances were then calculated.These data processing steps were performed in R version 4.0.5 [55], the R code is available in the GitHub repository https:// github.com/ rksal is/ fungi cides_ metab arcod ing and the Zenodo repository https:// doi.org/ 10. 5281/ zenodo.77973 89.For the PacBio long reads, samples were assembled with LAA version 2.4.2 [56].

Phylogenetic reconstruction
Phylogenetic relationships were reconstructed for a chosen set of hyphomycete taxa that covered the majority of the ASVs.For each individual taxon, the equivalent subtree was extracted from the 18S-based reference tree maintained by the creators of the ghost-tree pipeline (ghost_tree_100_qiime_ver8_dynamic_02.02.2019,Fouquier et al. [57]).Corresponding ASVs, long readcontigs and reference sequences from UNITE were aligned using MAFFT version 7.310 7.310 (maxiterate 1000-globalpair, Katoh and Standley [58]) and trimmed with trimAl v. 1.4.1 (-gt 0.1, Capella-Gutiérrez et al. [59]).Phylogenies were reconstructed with FastTree version 2.1.10[60] using the subtrees derived from the ghost-tree reference tree as constraints.
The phylogenetic reconstructions were, furthermore, used to assign genera names for the ASVs for which taxonomic identifications from DADA2 were inconclusive.In these cases, maximum parsimony hidden state prediction as implemented in the R package castor version 1.5.7 [61] was utilised to assign genera names to internal nodes.These labels were propagated to unassigned ASVs if they were distanced no further than 0.2 substitutions per site from the closest sequence of the corresponding genus.

Assignment of functional groups
Fungal ASVs were grouped to the lowest taxonomic level assigned (generally species or genus level) using phyloseq version 1.30.0.Subsequently, aquatic hyphomycetes and chytrids were assigned to functional groups based on expert knowledge, and the remaining taxa were assigned to the following fungal trophic modes using the tool FUNGuild version 1.1 [62]: pathotrophs, saprotrophs and symbiotrophs.Taxa that remained unassigned or assigned to mixed trophic modes were classified as "other".The Functional Annotation of Prokaryotic Taxa database version 1.2.6 (FAPROTAX; Louca et al. [63]) was used to assign functions to the bacterial ASVs, these include functions related to biogeochemical processes.

Statistical analysis
To compare morphological and DNA-based taxon identification approaches for fungi, the ASVs representing aquatic hyphomycete taxa were identified and grouped to the lowest taxonomic level assigned (generally species or genus).This resulted in a total of five community data sets (Additional file 1: Fig. S1): (1) aquatic hyphomycete taxa identified morphologically using conidia (sporulation rates were calculated from this data set, data missing from cycle 1), ( 2) aquatic hyphomycete species identified via ITS metabarcoding, i.e., unique species names extracted from the ASV data, (3) the entire fungal ASV data set identified via ITS sequencing (with and without successful taxonomic assignment), (4) the subset of aquatic hyphomycete ASVs (also those without species assignment) and ( 5) the entire bacterial ASV data set identified via 16S metabarcoding.In addition, data on total bacterial density and decomposed leaf mass were available.For the comparison of aquatic hyphomycetes identified based on morphology vs ITS sequencing, only cycles 2 and 3 were included in the manuscript, as morphological data from the first cycle was missing.A similar comparison including ITS sequence data from cycle 1 is included in Additional file 1.
Taxa richness (hereafter called species richness) was calculated for data sets 1 and 2 (aquatic hyphomycete species identified via morphology and metabarcoding).For data sets 3, 4 and 5, ASV richness was calculated (entire fungal data set, aquatic hyphomycetes subset, bacteria data set).Generalised Linear Models (GLMs) with Poisson distribution were used to determine the effect of fungicide exposure (categorical variable with 2 levels: control and treatment) and cycle (categorical variable with 3 levels: 1, 2, and 3; except for data set 1, with 2 levels: 2 and 3) and their interaction on the species and ASV richness of all data sets.Furthermore, linear models with log-transformed total bacterial density as well as decomposed leaf mass as response were fitted with the same predictors, using type II ANOVAs with F tests to evaluate statistical significance (defined as p < 0.05).The same explanatory variables were used in linear models for individual aquatic hyphomycetes that were identified with both morphological and sequencing methods to assess potential differences in their response between methods.
Potential changes in community (species or ASV) composition were analysed using Redundancy Discriminant Analyses (RDA).Before analysis, community data were Hellinger transformed to standardise the data and circumvent problems associated with the use of Euclidean distance for ecological data [64].We calculated the proportion of variance explained for the individual and joint explanatory variables (i.e., cycle and fungicide exposure), including their interaction, to evaluate the individual and shared fraction of variance explained [65].The statistical significance of the RDAs was checked using permutational type III ANOVAs.To determine which fungal and bacterial ASVs differed with fungicide treatment in each of the three cycles, differential abundance testing was performed using the R package "DESeq2" version 1.38.3 [66], which uses negative binomial models with Wald statistics.Log2-fold change values were calculated and Wald test P values were adjusted for multiple testing using the Benjamini-Hochberg procedure.
All statistical analyses were performed in R version 4.0.5 [55].Richness calculations and RDA analyses were performed using the R package vegan version 2.6.4 and the GLMs/LMs were run using the base R functions.The raw data for this study have been uploaded to the NCBI Sequence Read Archive (SRA) under BioProject PRJNA944599 (ASV sequences) and PRJNA991450 (raw PacBio long reads) and to GenBank (OR243761:OR243775; long consensus sequences).The raw data and R code used for the analyses are available in the GitHub repository https:// github.com/ rksal is/ fungi cides_ metab arcod ing and the Zenodo repository https:// doi.org/ 10. 5281/ zenodo.77973 89.

Results
Sequencing yielded a total of 137,491-330,202 ITS reads and 118,040-345,780 16S reads per sample, which were assigned to 1882 fungal ASVs and 5798 bacterial ASVs, respectively.Of the 1882 fungal ASVs, 271 belonged to aquatic hyphomycetes.The taxonomic assignment of fungal ASVs in our study was enhanced through the generation of long-read reference sequence data for locally occurring fungal key species as well as by applying a phylogenetic placement approach.Five additional ASVs could be assigned to species level by adding the long-read sequences of the isolated cultures to the UNITE database [33] to assign taxonomy.A further 43 ASVs could be assigned to genus level based on the constructed phylogenetic trees (Additional file 1: Table S2, Figs.S2-S6).

Comparison of community composition using morphological and DNA-based methods
Overall, 22 hyphomycete taxa were identified morphologically and eleven via ITS sequencing in cycles 2 and 3 (Fig. 1).Eight additional taxa were identified with metabarcoding during cycle 1 (Additional file 1: Fig. S7).Of these, six species were identified via both, morphological identification as well as metabarcoding.Many others were, however, closely related (belonging to the same genus).
Six morphologically identified taxa were only identified to genus level.For four, several ASVs of the same genera were identified with metabarcoding (Fig. 1), while for the remaining two no ASVs of the same genera were detected, although related congeneric species were present in the UNITE database.Ten species identified morphologically were not detected with metabarcoding.Of these, three were present in the sequence database (Alatospora pulchella, Lemonniera cornuta and Lemonniera terrestris), while the remaining seven were absent from the sequence database, although for three species (Stenocladiella neglecta, Fontanospora eccentrica and Cylindrocarpon aquaticum), related congeneric species were present in the database.
Five taxa (species or genera) were detected via metabarcoding but not identified morphologically (Fig. 1).Of these, four were groups of ASVs assigned to the genus level but with no species assigned and all had matching genera identified via morphology.However, for the remaining species identified via metabarcoding (Anguillospora crassa) no congeneric species was identified morphologically.

Species richness and community composition
The aquatic hyphomycete species richness identified morphologically (data set 1) was neither affected significantly by treatment nor by cycle (Fig. 2a, Additional file 1: Table S3), though missing data for cycle 1 may have masked a response.In contrast to this, aquatic hyphomycete species richness identified through ITS sequencing (data set 2) significantly decreased from cycle 1 to cycles 2 and 3 (p < 0.001), from an average of 13 aquatic hyphomycetes species in cycle 1 to approximately seven in cycles 2 and 3 (Fig. 2b, Additional file 1: Table S3).Similar to data set 1, it also showed no response to fungicide exposure.Both methods identified a similar species richness in cycles 2 and 3, with an average of seven and eight species, respectively (Fig. 2a, b).
ASV richness of aquatic hyphomycetes and all fungi identified via sequencing (data set 3) was significantly affected by cycle (aquatic hyphomycetes: p < 0.001, all fungi: p < 0.001, Additional file 1: Table S3), decreasing from cycle 1 to cycles 2 and 3 (Fig. 2c, d, Additional file 1: Table S3).ASV richness also showed a significant interaction between treatment and cycle (aquatic hyphomycetes: p = 0.026, all fungi: p = 0.007, Fig. 2c, d, Additional file 1: Table S3).The fungicide exposure led to a lower richness in ASVs after cycle 1 in comparison with the control.This was, however, reversed in the following cycles, leading to a lower richness in the control than the fungicide treatment by the third cycle (Fig. 2c).
Of the six most common species found with molecular and morphological identification (data set 4), none showed significant responses to the fungicide treatment (Additional file 1: Table S5, Fig. S8).However, one species (Tetrachaetum elegans) decreased over cycles in both approaches (morphological: p = 0.004, molecular: p = 0.002, Additional file 1: Fig. S8, Table S5).In only considering molecular identification, Alatospora acuminata Fig. 1 Comparison between the aquatic hyphomycete community composition (in cycles 2 and 3) identified morphologically and using ITS sequencing.Colour intensity indicates the relative sporulation rate (in conidia per g leaf dry mass) and relative read abundance, respectively (see legend).Numbers in parentheses (right of ITS read abundance) indicate the number of amplicon sequence variants (ASVs) assigned to the species.*Articulospora tetracladia (anamorph) recorded in the UNITE database as Hymenoscyphus tetracladius (teleomorph).$ Indicates species identified morphologically for which there was no reference sequence in the database.The same figure including the ITS sequencing results of cycle 1 can be found in Additional file 1: Fig. S7 decreased over cycles (p < 0.001), while Articulospora tetracladia and Tetracladium marchalianum increased in relative abundance with cycle (p < 0.001 and p = 0.005, respectively; Additional file 1: Fig. S8, Table S5).
At the intraspecific level, in contrast, some ASVs (30 of 1882 fungal ASVs) showed significant responses to the fungicide treatment (Additional file 1: Table S6), including those which did not exhibit effects at the species level.Overall, in cycle 1, seven fungal ASVs decreased in their relative abundance with treatment and two increased; in cycle 2, six decreased and two increased; and in cycle 3, 14 ASVs increased with treatment, while no ASVs decreased (Additional file 1: Table S6).The majority of ASVs which showed a significant response toward the fungicide treatment were aquatic hyphomycetes.Importantly, among the ASVs responding to the treatment but assigned to other fungal groups, two were classified as saprotrophs, one was classified as a pathotroph and two more are known to exhibit both trophic modes (Additional file 1: Table S6).

Functional responses
Leaf decomposition showed a temporal stress response with a significant interaction between cycle and fungicide treatment (F value = 5.9, p = 0.006).This was mainly driven by an initial reduction in leaf decomposition by fungicide exposure (Fig. 2e).This effect decreased during the subsequent cycles with leaf decomposition matching control levels by the third cycle (Fig. 2e).The percentage of aquatic hyphomycetes in the fungal community increased with both fungicide treatment and cycle (treatment: p = 0.042, cycle: p < 0.001, Additional file 1: Fig. S9).

Responses of the bacterial community Richness and community composition
Fungicide treatment and cycle showed a significant interaction for bacterial richness (data set 5; p < 0.001, Fig. 2f, Additional file 1: Table S3).The bacterial density significantly responded to cycle (F = 5.7, p = 0.011, Fig. 2 g), increasing in density from cycle 2 to cycle 3.RDAs revealed that cycle was the major driver of the bacterial community composition, explaining 48.7% of the variation (p < 0.001, Fig. 3e, Additional file 1: Table S4).
Few significant responses (twelve of 5798 bacterial ASVs) to the fungicide treatment were also seen at the ASV level for bacteria (Additional file 1: Table S6).In cycle 1, four ASVs decreased with treatment and four increased; in cycle 2, one decreased and one increased; and in cycle 3, two ASVs decreased with the fungicide treatment (Additional file 1: Table S6).

Functional responses
The most common bacterial functions inferred using FAPROTAX were chemoheterotrophy, aerobic chemoheterotrophy, chloroplasts, and intracellular parasites followed by nitrate reduction (Additional file 1: Fig. S10).Assigned functions showed a significant response to cycle (F value = 15.945,p < 0.001) but not fungicide treatment (Additional file 1: Fig. S10).

Differences in aquatic hyphomycete community composition identified through morphological and molecular methods (H1)
We compared the aquatic hyphomycete community composition assessed with morphological identification of conidia to the molecular method of ITS metabarcoding with the primary focus on cycles 2 and 3 for which both morphological and molecular data were available.We hypothesised that most taxa found morphologically would be also found with DNA metabarcoding.Contrary to this expectation, we found a higher number of species with morphological identification compared to DNA metabarcoding (22 and eleven species, respectively).Accordingly only 22% (six species, Fig. 1) of all of the aquatic hyphomycete species detected in cycles 2 and 3 (17% when including ITS results from cycle 1, Additional file 1: Fig. S7) were found using both identification methods (i.e., morphological and molecular).This is a lower ratio of species identified with both methods than in a previous study (53%) [67].For both identification methods, however, several hyphomycetes could only be identified to genus level (morphological: 6, molecular: 4) as species level could either not be assigned or taxa could not be identified further.The number of species behind the four taxa identified via metabarcoding is unclear as these four genera consisted of a total of 16 ASVs (Fig. 1).When only taxa assigned to species level are considered the percentage of the aquatic hyphomycetes detected by both identification methods increases to 35% (26% when including cycle 1).The species identified by both methods in the present study were among the most dominant ones in the community and many other closely related species of the same genera were identified (Fig. 1).For example, while Orbilia rosea (previously Anguillospora rosea) was identified morphologically but not via metabarcoding (this species is absent from the reference databases), molecular data did reveal the presence of two morphologically similar species (Anguillospora crassa and Anguillospora longissima).Similarly, Lemonniera cornuta and Lemonniera terrestris could be identified morphologically to species level (with more spores only being identified to genus level), while seven ASVs could only be assigned to the Lemonniera genus (Fig. 1).Overall, all five of the taxa detected in cycles 2 and 3 via metabarcoding had genus-level counterparts identified morphologically, while only four of eight taxa identified in cycle 1 were identified to genus level with morphological identification of the later cycles.
Metabarcoding is generally limited by the marker choice and completeness of the used reference databases.The taxonomic assignment of fungal ASVs in our study was enhanced through the generation of long-read reference sequence data for locally occurring fungal key species as well as by applying a phylogenetic placement approach.An additional five ASVs were assigned to species level using the long-read sequences of isolated cultures and 43 ASVs were assigned to genus level using phylogenetic placement (Additional file 1: Figs.S2-S6).This highlights the potential of improving metabarcoding using longer reads [39] and using phylogenetic placement methods.Furthermore, metabarcoding could be improved by developing the databases further (which requires morphological identification in combination with sequencing), which is evident in the present study as seven of the ten species only identified morphologically were absent from the sequence database.By improving the databases, we may be able to capture a substantially greater proportion of the aquatic hyphomycetes/fungal community.This will ultimately allow us to take advantage of already collected metabarcoding data by re-running the analytical pipeline with improved databases, which, furthermore, highlights the need to publish all raw data to enable such reanalysis.Metabarcoding can, furthermore, capture dormant or even dead organisms, as long as the DNA remains intact [68].This could lead to misinterpretation of results, especially when investigating stress responses and being interested in the potential future community composition (i.e., healthy organisms which are able to reproduce) [69].One clear gain of using metabarcoding is that it enables the differentiation of species which can morphologically not be identified unambiguously as well as the differentiation of individual sequence variants (ASVs) allowing a better community resolution and, therefore, a higher explanatory power (Fig. 3).The fact that we were able to identify low abundant species only using metabarcoding but not morphologically (Fig. 1), indicates that metabarcoding may also enable the detection of rare taxa.By contrast, morphological identification can be pushed to its limits as cryptic species, phenotypic plasticity or other issues prohibit an unambiguous identification [29,30].Furthermore, morphological identification based on conidiospores, i.e., asexual reproduction products of the anamorph stage of fungi, as used here, captures only active organisms with the ability to reproduce under the present conditions and in the respective life stage [70,71].Another challenge of matching morphological and molecular identification is rooted in the fact that hyphomycetes can have two reproductive states, asexual (anamorph) and sexual (teleomorph) [72].In the past, this led to taxonomic confusion with two species names sometimes given to different stages of the same fungus (example faced in the present study: Articulospora tetracladia (anamorph) identified via morphological identification was matched to Hymenoscyphus tetracladius (teleomorph) as recorded in the UNITE database) and a unified taxonomy has not been achieved for all species yet.Previous studies, furthermore, demonstrated that species easily distinguishable in morphological identification have only barely differing ITS sequences with pairwise sequence differences of, e.g., 1.1% (Heliscus submersus and Flagellospora penicillioides) or 1.3% (Lemonniera alabamensis and Lemonniera aquatica) [73,74].
Overall, we found fewer hyphomycete species with metabarcoding and the overlap between metabarcoding and morphological identification was low, contradicting our first hypothesis.Metabarcoding, however, allowed to expand the investigated species pool by enabling the identification of species beyond hyphomycetes with identifiable spores and thus providing a higher explanatory power for the response within the experiment (Fig. 3).This is in line with the results of Pesce et al. [75], who were only able to detect community alterations of fungicides using molecular methods (namely, Automated Ribosomal Intergenic Spacer Analysis) but not morphological identification.

Comparison of functional and structural response patterns of the fungal community (H2)
As reported by Schreiner et al. [41], leaf decomposition initially declined strongly in response to fungicide exposure but recovered during the following two cycles (Fig. 2e).This may be explained by an adaption of the microbial communities to fungicide stress, which ultimately allowed the communities to recover their functional performance.In the present study, we used DNA metabarcoding to detect potential community changes that would reflect the functional response, which may have been missed using morphological identification of aquatic hyphomycetes via spores.
In contrast to our second hypothesis, we found similar community response patterns with both molecular and morphological identification at the species level (Fig. 3c,  d).Both methods showed that aquatic hyphomycete species richness did not respond to fungicide exposure.The species richness of aquatic hyphomycetes identified through ITS sequencing significantly decreased with cycle, whereas this was not detected morphologically, most likely due to missing data (Fig. 2a, b).Indeed, a companion study found such a decrease via morphological identification in communities from other biogeographic regions (i.e., Sweden and Denmark, Schreiner et al. [41]).Both morphological and sequencing methods identified a similar number of aquatic hyphomycetes in cycles 2 and 3 (seven to eight species), lending further support to this hypothesis (Fig. 2a, b).
The response of fungal communities to the cycle irrespective of the method used can likely be explained by adaption to laboratory conditions, in line with previous studies using field communities in microcosms [76,77].Differences between laboratory and field conditions, for example, in temperature and flow (Additional file 1: Text S1, Fig. S11), likely altered sporulation rates of various species [70].The loss of taxa may be explained by species with low sporulation rates being outcompeted by other taxa during the colonisation of new substrate [78].
While both methods, when compared at the species level, identified cycle as the only significant driver of aquatic hyphomycete community composition (identified via RDA, Fig. 3c, d, Additional file 1: Table S4), a considerably higher proportion of variance was explained when using metabarcoding (53% explained variance via molecular identification, compared to 17% for morphological identification).The clearly higher ratio of explained variance emphasises the advantages of molecular identification in comparison with morphological identification, such as including non-sporulating hyphomycetes or unambiguous identification among others.Although less pronounced, a similar trend of a higher explained variance when using molecular methods was previously found for hyphomycete identification based on denaturing gradient gel electrophoresis (DGGE) [79].
The pattern of decreasing richness with cycle was also observed when looking at ASV richness of aquatic hyphomycetes and all fungi identified via molecular identification (Fig. 2c, d).ASV richness of aquatic hyphomycetes and all fungi was significantly affected by an interaction between treatment and cycle, in contrast to species richness of hyphomycetes (Fig. 2c, d; Additional file 1: Table S3).The fungicide exposure initially led to a lower ASV richness (cycle 1), this was, however, reversed in the following cycles, leading to lower ASV richness in the control treatment compared to the fungicide treatment by the third cycle.These trends in ASV richness match the pattern in leaf decomposition (Fig. 2e, see above).Previous studies showed that hyphomycete species identity has stronger explanatory power for decomposition than richness [12,80].
Furthermore, while no significant effects of the fungicide exposure were detected when regarding the fungal community composition (Fig. 3a-d) or individual aquatic hyphomycete species identified morphologically or molecularly (Additional file 1: Table S5), significant effects on individual ASVs were observed (Additional file 1: Table S6, details on example species Alatospora acuminata in Additional file 1: Text S2).The higher taxonomic resolution achieved using molecular identification enables the detection of different responses of individual ASVs of a single species since even within single species, various genotypes (ASV-level) can exhibit different tolerances toward fungicides.
Although the ASV and species richness of all fungi and aquatic hyphomycetes in particular, decreased with cycle (Fig. 2a-d), the contribution of aquatic hyphomycetes in the fungal community (as detected via metabarcoding) increased across the subsequent cycles to nearly 100% (Additional file 1: Fig. S9).This suggests that hyphomycetes, which are known to be important for decomposition, can despite known challenges of adapting reproduce under microcosm conditions (Additional file 1: Text S1, Fig. S11 [76,77]), which is not the case for other aquatic fungi, as suggested by the presented data.Current research, however, is limited to studying aquatic hyphomycetes, which prohibits comparisons to other studies.
Overall, despite observing functional impairment in decomposition, especially in cycle 1 (Fig. 2e), little fungal community turnover was found (Fig. 3, Additional file 1: Table S4).Alterations in the function, i.e., decomposition, due to fungicide exposure were observed in multiple field as well as laboratory studies (e.g., [81][82][83][84]).These changes in leaf decomposition could regularly be explained by alterations in the fungal or bacterial community [81,83,85,86], although the fungicide exposure was mostly higher than in the present study.While studies imply that functional and structural responses are coupled, our results and those of Fernández et al. [79] only found functional but no structural changes and others only observed structural but not functional changes in response to fungicide exposure [87,88].Such varying results could, however, be caused by differences in the species composition of the fungal communities in the different studies.This highlights that structure and function can currently not be related reliably and that understanding the underlying mechanisms of the functional changes in response to fungicides but also other stressors requires further studies analysing both diversity and function in more detail.

Comparison of functional and structural response patterns of the bacterial community (H3)
Similar to the fungal communities ("Comparison of functional and structural response patterns of the fungal community (H2)" section), changes in the bacterial community composition as well as the total bacterial density were driven only by cycle (Fig. 3e, Additional file 1: Table S4, Fig. 2g), partially contradicting our third hypothesis.This lack of effect of fungicide exposure on the bacterial community composition is consistent with two previous studies [81,88], but contradicts another [79].Furthermore, a previous study detected increasing bacterial density throughout a long-term experiment but no effect of a fungicide treatment, which is in line with our findings [89], while another study detected decreasing bacterial density with increasing fungicide exposure [85].The bacteria ASV richness was, however, affected by both, cycle and the fungicide exposure (Fig. 2f, Additional file 1: Table S3), indicating complex bacterial community changes.From a total of 5798 bacterial ASVs, only a small minority (twelve, Additional file 1: Table S6) showed significant responses to the fungicide exposure.The majority of these responses (eight out of twelve) were detected during cycle 1 (Additional file 1: Table S6), which also exhibited the strongest changes in leaf decomposition (Fig. 2e).However, since as many bacterial ASVs increased in read abundance as decreased (four each in cycle 1), a causal relationship of these observations is unlikely (Additional file 1: Table S6).
The investigated bacterial communities exhibited a wide range of ecological functional groups related to the carbon cycle, such as chemoheterotrophy, aerobic chemoheterotrophy, phototrophy, photoheterotrophy, hydrocarbon degradation, aromatic compound degradation, methanol oxidisation and methylotrophy (Additional file 1: Fig. S10).Changes in these functional groups were only influenced by cycle and not by fungicide exposure (Additional file 1: Fig. S10).This is in line with a previous review, which stated the tolerance of these bacterial functional groups toward fungicide exposure although in soil [90].The dominating functional groups (i.e., chemoheterotrophy and aerobic chemoheterotrophy), related to the use of organic material as main sources of carbon and energy, were less common (expressed as mean number read numbers) than in cycles 2 and 3 (Additional file 1: Fig. S10).An increase in the relative abundance of these functional groups over cycles, together with other less common functional groups (e.g., methylotrophy or methanol oxidation) could be related to the increasingly similar decomposition in the treatments over cycles (Fig. 2e).

Conclusion
Morphological and molecular identification showed similar responses of fungal communities to two factors in a microcosm experiment, fungicide exposure and cycle.This, together with the fact that a similar number of hyphomycete species were identified, suggests that both identification methods provide reliable data when investigating hyphomycetes associated with leaf decomposition in streams.Consequently, results of studies working with molecular methods can be compared with results of studies using traditional morphological identification including those dating back decades.Future studies could, however, benefit from a higher explanatory power of molecular identification, as community turnover can be studied in much greater detail than with traditional methods.This also includes different tolerances of ASVs within one species or cryptic species, which might be the cause of contrasting responses of single species observed in previous studies.Further studies focussing on stress responses in terms of abundance, biomass and function (i.e., decomposition) of single ASVs as well as communities with known compositions are required to understand the mechanisms of alterations of microbial communities and their function.
tree of the order Pleosporales (fungi), constructed from sequences from the UNITE database [33], long read sequences of fungal cultures (indicated as orange bars) and ITS metabarcoding sequences from the experimental samples (indicated as purple dots).The outer colour strip indicates the families the sequences are assigned to (see legend on left-hand site).Figure S4.Phylogenetic tree of the order Hypocreales (fungi), constructed from sequences from the UNITE database [33], long read sequences of fungal cultures (indicated as orange bars) and ITS metabarcoding sequences from the experimental samples (indicated as purple dots).The outer colour strip indicates the families the sequences are assigned to (see legend on left-hand site).Figure S5.Phylogenetic tree of the order Tremellales (fungi), constructed from sequences from the UNITE database [33] and ITS metabarcoding sequences from the experimental samples (indicated as purple dots).The outer colour strip indicates the families the sequences are assigned to (see legend on left-hand site).Figure S6.Phylogenetic tree of the family Psathyrellaceae (fungi), constructed from sequences from the UNITE database [33] and ITS metabarcoding sequences from the experimental samples (indicated as purple dots).The outer colour strip indicates the families the sequences are assigned to (see legend on left-hand site).Table S3.Results of the generalised linear models (GLM) testing the influence of explanatory variables on responses observed via the two identification methods and the five data sets generated.Bold p values indicate statistical significance.Data was missing for the morphologically identified hyphomycete community from cycle 1.LRT: likelihood-ratio test.Table S4.Explained variance and p values of the Redundancy Discriminant Analysis (RDA) variables with Hellinger transformed data, tested by type III ANOVAs.Bold p values indicate statistical significance.Data was missing for the morphologically identified hyphomycete community from cycle 1.Table S5.ANOVA (type II) results testing the Influence of explanatory variables on the aquatic hyphomycete taxa richness observed via both identification methods, morphological and molecular.Figure S7.Comparison between the aquatic hyphomycete community composition (over all cycles) identified morphologically based on morphologically and using ITS sequencing.Colour intensity indicates the relative sporulation rate (in conidia per g leaf dry mass) and relative read abundance, respectively (see legend).Numbers in parentheses (right of ITS read abundance) indicate the number of ASVs assigned to the species.*Articulospora tetracladia (anamorph) recorded in the UNITE database as Hymenoscyphus tetracladius (teleomorph).$ indicates species identified morphologically for which there was no reference sequence in the database.Figure S8.Responses of four common aquatic hyphomycete taxa (that were identified via both morphological identification and metabarcoding and showed at least one significant response) to the different cycles (numbers on top) and treatments (circle: control; triangle: fungicide treatment).Mean relative abundance measured via sequencing is shown in the four upper plots and mean sporulation (morphological identification, in conidia per g leaf dry mass) is shown in the four lower plots.Error bars represent 95% confidence intervals.P values are provided for significant results based on the two-way ANOVAs (type II).Data was missing for the morphologically identified hyphomycete community from cycle 1.Table S6.Differential analysis of fungal and bacterial amplicon sequence variants (ASVs) using DESeq2 [66].Showing the log-fold change values of ASVs significantly different in treatment and controls for cycles 1, 2 and 3. Positive log fold change values means that ASV normalised read abundance increased in treatment compared to control, and negative log fold change value that ASVs read abundance decreased.Fungi marked with * are saprotrophs, with † are pathotrophs and with § are mixed (exhibit both trophic modes).

Fig. 2
Fig. 2 Species richness for fungi and bacteria, decomposed leaf mass and bacterial density across cycles and treatments.Mean species richness for the aquatic hyphomycetes based on sporulation and morphological identification (a), the ITS-identified aquatic hyphomycete species (b), the ITS-identified aquatic hyphomycetes amplicon sequence variants (ASVs) (c), and the entire ITS data set (d), the mean decomposed leaf mass [%] normalised for temperature (degree day: dday) (e), mean bacterial (16S) ASV richness (f), and bacterial density [log bacteria per cm 2 ] (g) for the different cycles (numbers on top) and treatments (circle: control; triangle: fungicide treatment).All mean values are displayed with 95% confidence intervals.P values are provided for significant results based on the Generalised Linear Models and two-way ANOVAs (type II)

Fig. 3
Fig. 3 Community composition across cycles and treatments.Community composition of the entire ITS data set (a), ITS-identified aquatic hyphomycetes amplicon sequence variants (ASVs) (b), ITS-identified aquatic hyphomycete species (c), aquatic hyphomycetes based on sporulation and morphological identification (d) and bacteria based on 16S metabarcoding (e).Each point within the Redundancy Discriminant Analysis (RDA) represents the community of one replicate.Replicates of one treatment are connected through lines (for molecular methods n = 10 including two technical replicates, for morphological identification n = 5).Data were missing for the morphologically identified hyphomycete community from cycle 1.P values and explained variance are shown for significant RDA variables tested by type III ANOVAs

Figure S9 .
Percentage of the different fungal groups (aquatic hyphomycetes and chytrids) and trophic modes (pathotrophs, saprotrophs, symbiotrophs and other unclassified fungi) in the different cycles and treatments.P values are provided for significant results based on the two-way ANOVAs (type II).
Figure S10.Bubble plot showing the most common bacterial ecological functional groups.The size of the bubbles is referring to the mean read numbers per treatment and cycle, while the colour of the bubbles is referring to the treatment (control = blue, fungicide = yellow).Text S1.Differences field and lab conditions.Figure S11.Experimental setup of the microcosm experiment.Flow was created using magnetic stirrers.The magnetic stirrer bar was separated from the leaf bags using a stainless-steel mesh.Photo by Verena C. Schreiner.Text S2.Response on the ASV level.