A metatranscriptomic exploration of fungal and bacterial contributions to allochthonous leaf litter decomposition in the streambed

Experimental setup and sampling

The experiment was conducted at the greenhouse of the University of Duisburg Essen from November 11th to 30th where six Aquaflow mesocosms (Graupner et al., 2017) were used (Fig. S1). Each system had two large tanks which were connected by two sediment-filled channels (narrow and broad). Each system was filled with 365 L of water and 60 L of sediment from the Boye River (51°33′17.8″N 6°56′43.2″E). Water was pre-filtered with a nylon net bag (200 μm) to keep out metazoa and sediment was mixed with a concrete mixer to homogenize the material. After homogenization, 45 cylindrical sediment cores (Height: 9 cm and Diameter: 3 cm) were placed in the channels of each system in a single row. The remaining portion of the flume was then filled with sediment to a height of 10 cm. Alder (Alnus glutinosa) trees were purchased from Eggert Baumschule (Vaale, Germany) and cultivated in the botanical garden of the University of Duisburg-Essen. Senescent leaves were collected from alder trees in the autumn of 2021. Only leaves were taken that were easily detached because they had already developed a parting tissue, indicating that they were about to be naturally abscised. The collected leaves were subsequently air-dried for a week. Dry leaves were put in small bags of garden net (5 mm). Each bag contained 0.8 g of leaves (~2–3 leaves) without a petiole. The collected leaves were also used to cut leaf-discs (diameter 10 mm) that were placed in fine mesh bags (36 bags, three discs per bag). All leaf bags were incubated with mixed leaves from the Boye River bed before the experiment to ensure colonisation by fungi and bacteria. For incubation, 400 partially decomposed leaves from Boye River such as hazelnut, poplar, birch, and alder were mixed with the alder leaf bags with 180 L of water from the Boye River. Three bubble air pumps were placed for oxygen supply at the bottom of the tanks. Following a seven-day incubation, the leaf litter bags and discs were introduced into the narrow channels of three systems (treatment), while the remaining three systems were left without leaves (control). Each channel contained 90 bags of leaves and six bags of leaf discs, placed above the sediment cores. Two large bags (25 × 20 cm), each containing 25 leaves, were placed in large tanks for initial colonisation in each of the six systems from the start of acclimatisation. These bags were removed after the fourth day of acclimatisation. The water temperature in each channel was set to 15 °C with the help of water-refrigerator TC20 (TECO SRL, Ravenna, ITA). This setup was acclimatised for 10 days by exchanging the water between all six systems daily to obtain the same initial conditions. Following the acclimatisation phase, sediment was systematically sampled for 10 days, whereas the leaf discs were sampled on days 1, 5, and 10. For sediment sampling, four cores were extracted from each flume every day and then homogenized together in a beaker. Subsequently, 2 g of the homogenized sediment was aliquoted into five separate 2 ml Eppendorf tubes. A total of 60 sediment samples were collected. These samples were promptly stored in a cryo-shipper charged with liquid nitrogen for transportation to the laboratory of the University Duisburg-Essen. They were kept at −80 °C until the commencement of nucleotide extraction.

Nucleic acid extraction and sequencing

To extract RNA from sediment samples, 5 g of sample material was processed as described in Buchner & Wolany (2023). Briefly, samples underwent homogenisation and lysis using an 800 μl guanidine isothiocyanate (GITC) lysis buffer and phenol-chloroform, employing a bead-beating method. RNA precipitation and binding were carried out using silica spin columns. Following thorough washing, RNA was eluted with 100 μl elution buffer at room temperature. The obtained RNA was subsequently subjected to sequencing at CeGaT NovaSeq in Tübingen, Germany. Library preparation and rRNA removal were achieved using Illumina TruSeq Stranded RNA and Ribo-Zero kits. Samples were sequenced with CeGat NovaSeq 6000, 101 bp sequencing (2 × 101).

For amplicon sequencing, DNA was extracted from 500 mg of each sample for bacteria and fungi. A mixture of 0.1 and 0.5 mm diameter glass beads and 100 μl Proteinase K, 5 μl RNAse A, and 900 μl TNES were used for bead beating method with Mini-Bead-Beater 96 (Biospec Products, Bartlesville, OK, USA) for 2 min at 2,400 rpm (for buffer and reagents see materials in Buchner (2022a). Samples were incubated at 56 °C and bead beaten at 1,400 rpm for 20 min. PCR replicates were divided from lysates and the spin column protocol using a vacuum manifold was used for DNA extraction (Buchner, 2022b). DNA clean-up (Buchner, 2022c) was performed with 40 μl DNA input and 80 μl of clean-up solution. A two-step PCR approach was used for DNA amplification. For the amplification of bacterial DNA, 515f/806r primers (Apprill et al., 2015) were used with a reaction volume of 10 μl, incorporating 1 μl of DNA. Similarly, for fungal DNA, fITS7 (Ihrmark et al., 2012) and ITS4 (White et al., 1990) primers were used with a 10 μl reaction volume and 1 μl of DNA in the initial PCR (Multiplex PCR Plus Kit; Qiagen, Hilden, Germany). First PCR is followed by the cleanup (Buchner, 2022c) and then 2 μl used for the second PCR. The cycling condition for the first and second PCR is given in Table S1. Sample where noramlised to ~2 ng/μl a with bead-based normalisation (Buchner, 2022d). For pooled libraries, a Spin-column clean-up protocol was used which resulted in a final volume of 100 µL (Buchner, 2022b). Two PCR replicates were produced from each sample to ensure accuracy and reliability. These libraries were then subjected to paired-end sequencing for bacteria (2 × 250 bp) on Illumina NovaSeq and fungi (2 × 300 bp) on Illumina MiSeq V3 (CeGat Gmbh, Tübingen, Germany). Amplicon sequencing data were generated utilizing an established, automated pipeline developed by a large consortium handling numerous samples and projects. In contrast, metatranscriptome sequencing required phenol-chloroform-based extractions to produce sufficient quantities of RNA. This extraction method could not be incorporated into the automated pipeline. However, as both extraction protocols utilize SDS and bead-beating for cell lysis, we consider the resulting data to be comparable.

Bioinformatic analysis and statistics

Effect of allochthonous leaf litter on global expression patterns

Raw reads were quality-filtered, assembled, and annotated with functions and taxonomy as mentioned in the material and methods. PCA revealed a slight separation between the control (without leaves) and treatment (with leaves) channels as well as time points for both bacteria and fungi. However, this separation corresponded to a very low explained variance on the higher principal components for treatment (bacteria: PC 5, 3.54%; fungi: PC 6, 3.59%) (Figs. S2 and S3) and for time points (bacteria: PC 4, 4.04%; fungi: PC 4, 4.61%) (Fig. S4). Furthermore, a parametric multivariate statistical permutation test (PERMANOVA conducted using adonis2) indicated no significant differences (p-value ≥ 0.05) between the control groups and treatment groups for either bacteria or fungi. Additionally, there were no significant interactions between day and treatment effects (Table S3).

Effect on functional diversity

The richness and evenness of all KO counts were examined using Shannon diversity. Significant differences between treatment and control samples were observed on Day 4 (ANOVA: p-value = 0.042), Day 6 (ANOVA: p-value = 0.008), Day 9 (ANOVA: p-value = 0.047) for bacteria (Fig. 1A) and Day 4 (ANOVA: p-value = 0.03) for fungi (Fig. 1B). While the difference was not statistically significant, functional diversity index was mostly higher for treatment than for control at all time points, especially for fungi (Fig. 1). A closer look at these latter timepoints showed significantly higher expression of CPOM degradation enzyme in treatment group compared to the control. Notable examples include cellulase, xylanase, pectin lyase, and various glycosidases (Table S5).

Temporal dynamics of CPOM degrading enzymes

To analyse the temporal change of transcript abundances of enzymes, TCseq analysis was performed for all KOs, and clusters with ascending trends (upregulation over time) were used for further analysis (Fig. 2). For bacterial KOs, these clusters include a group of enzymes involved directly or indirectly in CPOM degradation such as 4-hydroxyphenyl acetate 3-monooxygenase involved in the degradation of aromatic compounds potentially derived from lignin breakdown, Formaldehyde ferredoxin oxidoreductase involved in the oxidation of aldehydes which can be part of lignin breakdown, cellobiosyl-(1->3)-beta-D-glucosidase; involved in the breakdown of cellulose, an important component of plant cell walls, and β-Glucosidase; responsible for breaking down glucosides into glucose, a crucial step in cellulose degradation (Table S7). In fungal KOs, clusters with an upregulation trend include enzymes such as glucuronoxylanase which is involved in the hydrolysis of xylans, a main component of hemicellulose found in plant cell walls and β-glucosidase, an enzyme involved in the hydrolysis of terminal, non-reducing β-D-glucosyl residues with the release of β-D-glucose. It is crucial in the degradation of cellulose into glucose, directly contributing to the degradation of plant matter (Table S7).

For certain genes, bacteria and fungi showed similar mean gene expression ratios over time between treatment and control. Selected KOs specific for CPOM degradation (Table S4) are shown in Fig. 3.

Differential expression analysis

Differential expression analysis with DESeq2 design ‘treatment*day’ showed that in bacterial KOs, there were 230 significant KEGG orthologs for the treatment effect (addition of leaf litter), 85 KOs for the time effect, and 102 for the interaction effect of treatment with time (Table S8). Figure 4 illustrates the significant differential expression of enzymes associated with CPOM degradation. Notably, several bacterial enzymes are highlighted, including 3-hydroxybenzoate 6-monooxygenase, which breaks down aromatic compounds from lignin, and sphingolipid C4-monooxygenase, involved in sphingolipid metabolism. Additionally, beta-galactosidase hydrolyzes beta-galactosides into monosaccharides, while mannan endo-1,4-beta-mannosidase contributes to the degradation of the plant cell wall component mannan. Peptidases are also present among other enzymes (Table S8). In addition, CAZyme glycoside hydrolase families such as GH20 and GH26 annotated with significant KOs also showed upregulation with treatment effect (Table S8).

In the case of fungi, 618 KEGG orthologs were significant for the treatment effect, 77 KOs for the time effect, and 96 KOs for the interaction effect of treatment over time (Table S8). Here, the significant KOs annotated to multiple enzymes were involved directly or indirectly in CPOM degradation such as 3-hydroxybenzoate 6-monooxygenase for aromatic compound degradation which are components of lignin, L-xylulose reductase Involved in pentose catabolism, important for the breakdown of hemicelluloses like xylan present in leaf litter, glucarate O-hydroxycinnamoyltransferase catalyzes the transfer of hydroxycinnamoyl groups to glucarate, related to the metabolism of plant-derived acids, Beta-galactosidase hydrolyzes beta-galactosides into monosaccharides, Alpha-mannosidase catalyzes the hydrolysis of alpha-mannosides to release mannose which is important for degrading mannans and glycoproteins in leaf litter are among other enzymes which are involved in carbohydrate metabolism, lipid metabolism, and protein degradation are listed in Table S8. In addition, significant KOs were also annotated to the CAZyme family Glycoside hydrolase: GH26 and GH35. Glycoside hydrolases like GH26 and GH35 are known for their involvement in the degradation of complex polysaccharides into simpler sugars, facilitating carbon release from organic matter. GH26 is typically associated with the hydrolysis of hemicellulose components, whereas GH35 is involved in breaking down Beta-galactosides. We conducted a KEGG pathway enrichment analysis to gain further insights into the pathways of upregulated differentially expressed genes. In the context of degradation, significant KOs (p-value < 0.05) in both bacteria and fungi are enriched in crucial pathways including starch and sucrose metabolism, glycan degradation, alanine-aspartate and glutamate metabolism, and other pathways involved in breaking down of complex organic compounds (Table S9).

Community composition based on amplicon sequencing

Amplicon analysis using the Natrix2 workflow and filtering (count data Table S10) generated 5295 OTUs from total 200 million raw sequences for bacteria and 108 OTUs from total 9 million raw sequences for fungi. To see the effect of treatment in the community composition ordination was performed with PCA showing minor separation between control and treatment for bacteria (principal components accounted for axis 1 12.86% and axis 2 7.71% of the total variance) and fungi (principal components accounted for axis 1 22.23% and axis 2 7% of the total variance) (Fig. S5). PC analysis also showed separation of samples with time points at lower axes for bacteria (PC6: 3.71%) and fungi (PC6: 4.33%) (Fig. S6). Observed differences were validated using a PERMANOVA test, which showed that fungi had significant changes in composition between control and treatment while bacteria showed significant differences for condition and day (Table S11). In addition, we observed a significant change in diversity between control and treatment samples for fungi (Fig. S7). A closer look at microbial composition revealed that bacterial communities Pseudomonadota and Bacteriodota were the most abundant phyla (Fig. 5A) but differential analysis of OTUs abundance on phylum-level showed no significant change with treatment. Fungal communities were predominantly composed of phyla such as Basidiomycota and Ascomycota (Fig. 5B) but none were related to aquatic hyphomycetes (see Supplemental Information) on the genus level (Franco-Duarte et al., 2022) which can be due to primer bias and usage of more universal primers. Basidiomycota and Ascomycota abundances also show a significantly higher abundance within treatment compared to control (padj < 0.05 & log2foldchange > 2) (Table S12).

Cooperative functional relationship between bacteria and fungi

Shared temporal expressions of beta-galactosidase, differential expressions of mannan endo14-beta-mannosidase, and KEGG pathways involved in breakdown of complex carbohydrate and aromatic compounds by both bacteria and fungi suggest potential synergism for organic matter degradation in sediment which is in line with the study by Wohl & McArthur (2001) in water of freshwater streams. In comparison to bacteria, fungi show more differentially expressed enzymes related to the breakdown of complex carbohydrates with the increase of organic matter which correlates to more fungal activity over time. In comparison to bacteria, fungi show more differentially expressed enzymes related to the breakdown of complex carbohydrates with the increase of organic matter which correlates to more fungal activity over time due to higher expression or higher fungal abundance (Gulis & Suberkropp, 2003; Romaní et al., 2006). The observed higher fungal activity aligns with their recognized primary role in the early stages of decomposition, where they initiate the breakdown of complex organic matter, making it more accessible for subsequent microbial processing.

Consistent bacteria amidst diverse fungi

The analysis of community composition using 16S rRNA gene sequences showed a stable bacterial composition with phyla such as Pseudomonadota, and Bacteriodota, being the most abundant. These are known to play a crucial role in the degradation of coarse particulate organic matter (CPOM) as mentioned in previous studies (Newton et al., 2011; Besemer et al., 2012; Purahong et al., 2016; Hayer et al., 2016). Similarly, the analysis using the ITS region reveals the presence of fungal phyla such as Basidiomycota, and Ascomycota which also showed significantly higher abundance in differential analysis. These fungal groups are known for their relevance in decomposition processes, as documented in studies by Meng et al. (2024), Lundell, Makela & Hilden (2010), and Kuramae et al. (2013). However, none of the phyla detected belong to a known aquatic hyphomycetes species (Franco-Duarte et al., 2022). The absence of this important group of freshwater decomposers can be experimentally explained by the physico-chemical conditions in the sediments, such as oxygen depletion and low turbulence (Chauvet et al., 2016; Cornut et al., 2014), usage of universal primer due to varying sequence lengths in hyphomycete’s ITS regions or potentially missing references of undetermined hyphomycetes in the Unite database.

However, we did not see significant changes for all days (PERMANOVA test) but we did observe significant changes in specific enzymes (inferred from DESeq2) and diversity and richness of CPOM-related enzymes for certain days (Fig. 1). Limitations of our studies could be a shorter duration of the study compared to other traditional experiments, as noted in studies like Hättenschwiler, Tiunov & Scheu (2005) and a long period of decomposition process by microbes (90 days) as mentioned in Petersen & Cummins (1974). Furthermore, the most important part of the CPOM decomposition occurs directly on the leaf litter (Cornut et al., 2014) which could not be comparatively assessed in flumes without leaves. Despite the different significance levels for various enzymes, the overall expression patterns of CPOM related KOs in both fungi and bacteria display a similar trend over time (Fig. 3). This similarity may indicate potential synergism which is also mentioned by Bengtsson (1992) and functional redundancy.