Biogeochemical Regimes in Shallow Aquifers Reflect the Metabolic Coupling of the Elements Nitrogen, Sulfur, and Carbon

Principal-component analyses of hydrochemical data and available molecular data (Fig. 1A), meaning metagenome taxonomic and functional profiles (Fig. 1B to D), revealed a grouping according to previously identified hydrochemical clusters (13) or monitoring wells. The ordination of analyzed hydrochemical data (Fig. 1B and S2) showed a closer linkage of cluster 2 (H3-2 and H4-1) and cluster 5 (H5-2) monitoring wells, compared to the distinct placement of cluster 3 (H5-1) and 4 (H4-2 and H4-3) wells. The use of only redox-relevant parameters (Fig. S3) revealed a separation according to clusters as well. A comparison of all three ordinations indicated a clear separation of wells H3-2 and H4-1 according to site in the case of the metagenome profiles, although these wells belong to the same hydrochemical cluster. Well H4-1 (cluster 2) was more distinct based on the taxonomic than functional metagenome profile. The ordination based on the metagenome taxonomic profile set wells H5-2 and H4-1 distantly apart. The functional profile (Fig. 1D) revealed, similar to the hydrochemical data (Fig. 1B), a distant grouping of cluster 4 (H4-2 and H4-3) wells. Quantitative PCR targeting the 16S rRNA gene as a proxy for cell numbers revealed copy numbers ranging from 5.33 × 10⁶ copies liter⁻¹ (PNK66, H4-3) to 2.63 × 10⁸ copies liter⁻¹ (PNK66, H5-3). Most wells showed copy numbers in the region of 10⁷ copies liter⁻¹ (Table S1).

The identification of the main drivers responsible for the observed clustering in underlying taxonomic and functional profiles showed that the clustering of the wells based on taxonomy is driven by the differential abundance of Nitrospirae (contribution PC1, 84%) and, to a lesser extent, Planctomycetes (contribution PC2, 75%) and Proteobacteria (contribution PC2, 19%) (Fig. S4). From a functional perspective, the clustering of wells was mostly governed by gene functions linked to assimilatory nitrate reduction (contribution PC2, K00370 narG, >40%), nitrification (contribution PC2, K10944 to K10946 pmo-amoABC and K10535 hao, >5%), and anammox (contribution PC2, K20932 to K20934 hzs, >50%). Sulfate reduction-related genes (K00958 sat and K00394 aprA) were identified to play a less prominent role (Fig. S5).

The analysis of metagenome taxonomic profiles was complemented by a corresponding analysis of generated metatranscriptome data sets. These should be more responsive and thus eventually more resolving, taking into account that they cover not only those microbial community members that are merely present but also those that are actively partaking in community gene expression. Metatranscriptome taxonomic and functional profiles revealed a distinct clustering of the cluster 5 and 4 wells, respectively (Fig. S6).

Taxonomic profiling based on metagenome (taxonomic assignment rate, 43% to 58%) and metatranscriptome (31% to 56%) data identified a few dominating taxa in both cases (Fig. 2). Abundant phylum-level groups based on metagenome taxonomic profiles included Proteobacteria (32% to 50%), Nitrospirae (4% to 33%), Actinobacteria (5% to 10%), “Candidatus Omnitrophica” (3% to 9%), “Candidatus Patescibacteria” (including “Ca. Parcubacteria” and “Ca. Microgenomates”) (0.7% to 6.5%), and Planctomycetes (4% to 20%) (Fig. 2A). The comparable broad ranges of relative abundance are reflected by the peak abundances of individual phyla at particular wells. The two candidate phyla “Ca. Omnitrophica” and “Ca. Patescibacteria” showed highest abundances at H3-2, while Nitrospirae and Planctomycetes appeared to be overrepresented at sites H4-1 and H5-2, respectively. Only a minor proportion of sequences were assigned to archaea. Euryarchaeota showed the highest relative abundances (1.13 and 1.26%) at wells H3-2 and H4-3. Correlation analyses between selected hydrochemical parameters and metagenome taxonomic profiles showed strong correlations between oxygen availability and the abundances of Nitrospirae (correlation, 0.88; P < 0.001) and Bacteroidetes (correlation, −0.82; P < 0.001). Planctomycetes abundances were found to be linked to ammonium availability (correlation, 0.60; P < 0.05) (Table S2).

Metatranscriptome taxonomic profiles revealed high abundances of additional phylum-level groups (Fig. 2B), including Bacteroidetes and, at well H4-1, Thaumarchaeota. Thaumarchaeota showed clear differences in relative abundance in a comparison of data from different sampling campaigns. Proteobacteria showed increased abundances in metatranscriptome data sets compared to metagenome data sets. “Ca. Patescibacteria,” “Ca. Omnitrophica,” and “Ca. Rokubacteria” were much less abundant in metatranscriptome data sets (Fig. 2B). Calculation of the ratios between metatranscriptome and metagenome data sets (Fig. 2C) revealed increases in relative abundances for Proteobacteria and Bacteroidetes by factors of 1.1 to 1.7 and 1.1 to 6.6, respectively. Abundances of “Ca. Patescibacteria” and Nitrospirae were reduced by factors of up to 5 and 40, respectively.

A more resolved taxonomic profiling on the order and family levels showed that a small set of taxonomic groups is responsible for the high abundances of the aforementioned phyla (Fig. S7). Abundant groups within the Bacteroidetes included Flavobacteriia and Chitinophaga. Nitrospira sequences were mostly affiliated with the family Nitrospiraceae. Various groups within the Proteobacteria showed peak abundances at particular wells (Fig. S7). Rhizobiales were abundant in metatranscriptome data sets, with values ranging from 2.8% to 7.1% at wells H4-3 and H5-1 (Fig. S7B). High relative abundances were seen for Sphingomonadales in metatranscriptome data sets (wells H4-3 and H5-1, 1% to 12.4%). Burkholderiales were overall highly abundant, featuring increased relative abundances at cluster 4 wells (H4-2 and H4-3), and Nitrosomonadaceae abundances were highest at well H4-1. Increased abundances of Desulfuromonadales and Acidiferrobacteraceae were observed for well H5-2 (1.1% to 2.2% and 1.8% to 5.8%, respectively). Thaumarchaeota sequences were to a large extent affiliated with Nitrosopumilaceae (Fig. S7B).

Functional metagenome and metatranscriptome profiles (assignment rates, 18% to 22%) revealed comparable high abundances of sequences being affiliated with nitrogen cycling (Fig. 3A and 4A). Screening of the defined functional modules for nitrogen cycling and denitrification showed high relative abundances of sequences linked to narG and narH (assignment rates, 0.1% to 0.2%) in metagenome data sets. In metatranscriptome data sets, the relative abundances for these functions were 5 times higher. Comparing the two aquifer assemblages with respect to nirK and nirS revealed increased numbers of metagenome sequences nirK in case of clusters 1 to 3 (assignment rates, 0.01% to 0.04%), while the same was true for nirS and clusters 4 to 5. nirK was the dominantly expressed type of nitrite reductase. The highest expression values were seen for wells H5-1 and H5-2 (Fig. 4A). Metagenome sequences could be assigned to all essential genes for complete and incomplete denitrification over the range of sampled monitoring wells, indicating the encoded potential for this metabolic process along the Hainich groundwater monitoring transect. In the context of nitrification, peak relative abundances for sequences assigned to amoABC and hao were detected for H4-1 and, to a lesser degree, H5-1. Increased relative abundances of hao were also seen for H5-2. Similar patterns were seen in metatranscriptome data. Sequences linked to anammox gene functions (hdh and hzs) were detected throughout the transect, with strongly elevated abundances at site H5-2. The same observation was made on the metatranscriptome level. Metatranscriptome data revealed nitrification to be the overall dominant process in terms of expressed gene functions, with relative abundances of >1% for amoABC and hao at site H4-1 (Fig. 4B). The proportions of sequences assigned to nitrogen fixation were comparatively low in the metagenome and metatranscriptome data sets (0.004 to 0.018%). Gene functions linked to sulfur metabolism showed inconsistent relative abundance patterns along the transect. Functions associated with assimilatory sulfate reduction (encoded by cysCDHIJNNC) and sulfur oxidation-related gene functions showed no obvious site-dependent abundance pattern. Increased relative abundances of up to 0.01% were detectable for gene functions linked to dissimilatory sulfate reduction (sat, aprAB, and dsrAB) in wells H3-2, H4-2, H5-1, and H5-2. In comparison to nitrogen cycling gene functions, sulfur metabolism was underrepresented at the metatranscriptome level.

In order to identify direct links between taxonomy and function, sub-data sets were extracted from metagenome and metatranscriptome data based on determined functional annotations and were subsequently taxonomically assigned (Fig. 3B and 4B). On the metagenome level, denitrification was found to be strongly dominated by Nitrospiraceae. Nitrosomonadaceae, Planctomycetaceae, Ectothiorhodospiraceae, Caldilineaceae, and Thermaceae contributed to the denitrification sequence pool to a much lower degree. Nitrospiraceae also dominated metatranscriptome data sets, but Nitrosomonadaceae and Nitrosopumilaceae showed notable abundances as well (Fig. 4B). Nitrosopumilaceae (3.78%) and Nitrospiraceae (7.98%) had a sizable share of nitrification-related transcripts, while Nitrosomonadaceae (84.10%) were highly overrepresented. Sulfate reduction-related sequences were primarily linked to Gallionellaceae, Acidiferrobacteraceae, and Peptococcaceae. Desulfovibrionaceae-related sequences made up a significant proportion in metagenome data, while we did not identify them as contributing to community gene expression. No obvious link was seen between Bacteroidetes and any of the assessed biogeochemical functions, although this phylum showed high abundances in metatranscriptome data sets (Fig. 2).

We determined relative abundances in metagenome and metatranscriptome data for gene functions linked to common carbon fixation pathways (Fig. 3 and 4). We observed relative abundances of 0.002% to 0.009% and 0.02% to 0.04% for rbcS and rbcL, respectively, both representative for the Calvin-Benson cycle on a metagenome level. Sequences assigned to the 3-hydroxypropionate-4-hydroxybutyrate cycle were hardly detected. The cluster 4 and 5 monitoring wells featured increased relative abundances of sequences linked to the Wood-Ljungdahl pathway. Selected gene functions of the reductive tricarboxylic acid (rTCA) cycle revealed increased relative abundances in the case of aclA and aclB at the cluster 2 and 3 and frdA at the cluster 4 and 5 monitoring wells. Metatranscriptome data sets showed increased abundances for gene functions linked to the Calvin-Benson cycle at site H3-2 (up to 0.92%) and to the Wood-Ljungdahl pathway at site H5-2 (up to 0.52%). The small pool of sequences associated with the 3-hydroxypropionate-4-hydroxybutyrate cycle was found to be almost exclusively made up by sequences linked to the Thaumarchaeota. Sequences affiliated with the Wood-Ljungdahl pathway were primarily encoded and expressed by Planctomycetes at sites where anammox-related gene functions were overrepresented in metagenome and metatranscriptome data sets. At sites with increased relative abundances of gene functions linked to carbon fixation, for instance, H5-2, it was possible to detect all relevant gene functions for respective pathways, for example, the Wood-Ljungdahl pathway.

Independent from assessing either metagenome or metatranscriptome data sets, between 1.0 and 1.7% of the queried sequences were assigned to carbohydrate-active enzymes (CAZymes) (Fig. 5A). The distribution of assigned sequences among defined functional modules (CAZymes linked to cellulose, chitin, xylan, and pectin turnover) revealed that the majority of sequences were assigned to cellulase (EC 3.2.1.4; metagenome, 7.5% to 12.3%; metatranscriptome, 8.1% to 11.4%), chitinase (EC 3.2.1.14; metagenome, 7.6% to 9.6%; metatranscriptome, 6.3% to 12.4%), and endoxylanase (EC 3.2.1.8; metagenome, 6.0% to 9.5%; metatranscriptome, 7.5% to 11.5%) functions. The relative abundances of all CAZymes showed only minor fluctuations independent of the sampling site and well. This was not only observed for metagenome but also metatranscriptome data sets. Taxonomic groups containing and actively expressing genes with functions linked to cellulose, chitin, and xylan turnover were determined to identify organisms potentially involved in complex polysaccharide breakdown in the Hainich groundwater (Fig. 5B). In all cases, taxonomic profiles revealed high relative abundances of Proteobacteria, Nitrospirae, and Bacteroidetes. High abundances of Nitrospirae were primarily seen for well H4-1 on the metagenome level. Bacteroidetes were in particular abundant in the upper aquifer assemblage. Planctomycetes featured increased relative abundance at H5-2. Lesser-represented phyla comprised Actinobacteria, Acidobacteria, and Verrucomicrobia. In metatranscriptome data sets, these phyla were hardly detected at all. Metatranscriptome data sets showed, in general, higher abundances of Bacteroidetes, while those of Nitrospirae were lower.

In a comparison of the clustering of monitoring wells according to hydrochemistry and molecular profiles, our data show that microbiome composition allows these monitoring wells to be distinguished more strongly than hydrochemical parameters. Multiparameter analysis with and without redox-relevant parameters revealed no difference in the grouping of the monitoring wells (13). Of the factors considered for hydrochemical clustering, those easily affected by redox processes are of high relevance from a microbial perspective, since they are commonly used as electron donors and acceptors for energy metabolism. However, multiparameter analysis with and without redox-relevant parameters revealed no difference in the grouping of the monitoring wells (13). The differences between hydrochemistry and molecular profiles in discriminating monitoring wells are thus presumably a result of the spatially restricted stimulation of microbial activity by the limited availability of (in)organic energy sources.

Multiparameter clustering grouped wells H3-2 and H4-1 closely together (13), while these wells were clearly separated based on metagenome taxonomic and functional profiles (Fig. 1B and C). Given that H3-2 and H4-1 are situated within distinct aquifer assemblages that primarily differ in terms of oxygenation, differences based on metagenome profiles can be attributed to the oxygen requirements of present microorganisms and processes carried out by them. This is true, for instance, when looking at the high relative abundances of Nitrospirae on the metagenome level (Fig. 2A) and Thaumarchaeota on the metatranscriptome level (Fig. 2B), both groups that are supposedly linked to nitrification, which requires oxic conditions and NH₄⁺ as an electron donor. In concordance with high relative abundances of nitrifying groups in the present and active pool of the microbiome, we observed high relative abundances for metagenome/metatranscriptome sequences affiliated with nitrification. The elevated levels of NO₃⁻ for H4-1 (13) are thus presumably a result of high nitrification activity. Previous work measured nitrification rates of 0.48 ± 0.09 and 0.64 ± 0.39 nmol NO_x liter⁻¹ h⁻¹ at H4-1 (14). In the case of H3-2, increased NO₃⁻ concentrations (13) are likely a consequence of vertically transported NO₃⁻ due to its high mobility or surface-derived NH₄⁺ that became oxidized along the groundwater flow path. Although HTU is suboxic to anoxic, the low oxygen concentrations at H3-2 apparently support the activity of common nitrifiers (19).

A comparison of the clustering based on metagenome and metatranscriptome profiles revealed rather different patterns (Fig. 1 and S3). The limited discrimination of monitoring wells on the metatranscriptome level was presumably due to an overrepresentation of Planctomycetes on taxonomic and anammox-related and nitrification-related gene functions on a functional level. A direct comparison of metagenome and metatranscriptome data sets was further impaired by data sets originating partially from different sampling campaigns and by differences in replication (Fig. 1A).

Observed differences between metagenomic and metatranscriptomic taxonomic profiles underline that the mere presence of organisms is a poor proxy for their activity in terms of gene expression. Common soil dwellers, such as Acidobacteria, Verrucomicrobia, Chloroflexi, and, by trend, Actinobacteria, were underrepresented in metatranscriptome data. Eventually relocated from their inherent habitat by vertical transport, respective groups can be considered surface echoes, limited in gene expression after being washed out and probably adapted to soil rather than aquatic environments. Groups, such as “Ca. Patescibacteria,” or “Ca. Omnitrophica,” which made up substantial proportions of the Hainich groundwater microbiome, also only marginally contributed to community gene expression (Fig. 2). The commonly assumed small sizes of “Ca. Patescibacteria” (9, 20, 21) and other microorganisms, such as Parvarchaeota (22–25), suggest that the intracellular space for ribosomes and thus, protein biogenesis, is limited, as the cell volume is mostly occupied by the genome. A limited protein biosynthesis machinery eventually requires a smaller mRNA transcript pool, which could in part explain the observed dilution of “Ca. Patescibacteria” abundance in a comparison of metagenome and metatranscriptome data. Other reasons for low abundances on the metatranscriptome level include metagenome profiles being partially derived from dormant and dead cells or extracellular DNA (eDNA). Although eDNA is an excellent nutritional source, since it provides carbon, nitrogen, and phosphorus, it can show prolonged half-life times in the environment if it gets stabilized on colloid or mineral surfaces (26–29).

A few taxonomic groups featured pronounced overrepresentations in the metatranscriptome. These comprised Thaumarchaeota, Bacteroidetes, and Planctomycetes. A related question is why their high community gene expression levels were not reflected in microbiome composition. It is tempting to speculate that environmental conditions reflected by available carbon (inorganic versus complex organic carbon) or electron donors (NH₄⁺) triggered their gene expression and thus, metabolic activity, while growth was limited by the comparable low energy output. It is a common misperception that high relative abundance equals high metabolic activity (30). Pester and colleagues (31) could show that low-abundance sulfate reducers (<0.01%) contribute significantly to sulfate reduction in a minerotrophic fen.

The Hainich groundwater microbiome is driven by nitrogen cycling. Using metagenomics and metatranscriptomes in a complementary way, we identified wells H4-1 and H5-2 as hot spots for nitrification and anammox, respectively (Fig. 6). These findings were in line with previous efforts addressing the functional diversity of the Hainich groundwater microbiome by biomarker analysis using PLFAs, amplicon sequencing, and metaproteomics (15, 17). Hainich groundwater microbiome PLFA profiles enabled the identification of three distinct PLFA populations affiliated with oxic/suboxic (H3-1, H3-2, H4-1, and H5-1), anoxic iron-rich (H4-2 and H4-3), and anoxic NH₄⁺-rich (H3-2, H5-2, and H5-3) groundwaters. This partitioning was reminiscent of the clusters identified by multiparameter analysis based on hydrochemistry (13). Oxic/suboxic groundwaters characterized by strong specific PLFA signals for Nitrospirae (17) were also found to be dominant on the metagenome level (Fig. 2). Besides, anoxic NH₄⁺-rich groundwaters with pronounced ladderance signals (17) showed a high abundance of anammox bacteria, based on our metagenome/metatranscriptome analyses (Fig. 2 and 4). A strong involvement of “Candidatus Brocadiales” in nitrogen cycling for well H5-2 was recently shown with available metaproteome data sets (15). Recent work determined anammox rates in the HTU to range from 3.5 to 4.7 nmol N₂ liter⁻¹ day⁻¹ (16).

High abundances of metagenome and metatranscriptome sequences linked to nitrification and anammox coincided with increased abundances of sequences linked to carbon fixation (Fig. 3 and 4). Nitrospiraceae, known to utilize the reductive tricarboxylic acid cycle for carbon fixation (32, 33), were the dominating organisms linked to nitrification on the metagenome level (Fig. 3), and we also saw elevated abundances of sequences linked to this carbon fixation pathway. The same pattern was seen for Nitrosomonadaceae. Their dominance for nitrification in metatranscriptome data coincided with elevated sequence abundances of Calvin-Benson cycle gene functions (Fig. 4). The increased relative abundance of Thaumarchaeota for nitrification on the metatranscriptome level was a consequence of this group participating in nitrification as well. Their mode of CO₂ fixation is a modified 3-hydroxypropionate/4-hydroxybutyrate cycle (34). Although gene functions linked to this pathway were found to be severely underrepresented, the expression of respective gene functions was exclusively linked to Thaumarchaeota at well H4-1. Anammox bacteria use the Wood-Ljungdahl pathway for inorganic carbon assimilation (35, 36), and their abundances at well H5-2 were matched by increased relative abundances for gene functions of this pathway (Fig. 3 and 4). The hot spot character of wells H4-1 and H5-2 for nitrification and anammox (Fig. 6) was a clear sign that chemolithoautotrophy in both aquifer assemblages is modulated by the overall limited availability of inorganic electron donors and acceptors.

Particular gene functions linked to nitrification and denitrification showed unexpected gene expression patterns and taxonomic affiliations. That was the case for hao at well H5-2 and for narG and narH at well H4-1. Hydroxylamine oxidoreductase (Hao) enzymes are generally linked to nitrification, where they catalyze the second step from hydroxylamine to nitric oxide/nitrite (37). Their increased expression at well H5-2 indicated a link with anammox. In fact, hao and hdh (hydrazine dehydrogenase gene) are considered homologs (38–40), which explains the unexpected increase in gene expression for this function. Our definition of functional modules (Fig. 3 and 4) does not include nxrAB in the case of nitrification; nxrAB encode nitrite oxidoreductase, catalyzing the oxidation of nitrite to nitrate. This is because of the joint grouping of nxrAB and narGH in two common gene functions (K00370 narG and nxrA and K00371 narH and nxrB) by KEGG. In fact, narGH and nxrAB are considered homologs due to high similarities at the amino acid sequence level (41). As a result, the increased expression of narGH at well H4-1 is presumably an increased expression of nxrAB and might be linked to nitrification (Fig. 4).

Metagenomics and metatranscriptomics allowed us to relate gene functions of different biogeochemical cycles to each other, thus bypassing a major constraint of functional marker-based analysis. Although nitrogen cycling-related aspects were highly overrepresented (Fig. 6), we observed an increase in sequences linked to dissimilatory sulfate reduction (sat, aprA, and dsrAB) at well H5-2 in metatranscriptome data, which was suggested by earlier PLFA profiles (17). Surprisingly, a high proportion of these sequences were annotated to Acidiferrobacteraceae (Fig. 4), a family recently described and only containing sulfur oxidizers (42). Screening of the Acidiferrobacteraceae genomes available from NCBI RefSeq revealed the presence of all genes essential for thiosulfate oxidation via the SOX system, as well as all genes needed for dissimilatory sulfate reduction and sulfide oxidation (sat, aprAB, and dsrAB) in Sulfurifustis variabilis (43). Although Acidiferrobacteraceae use the Calvin-Benson cycle for carbon fixation (43), we did not find Acidiferrobacteraceae-related sequences in our data. The anoxic conditions at well H5-2 rule out sulfur oxidation by oxygen reduction. Considering the presence of nirS and nirBD genes in S. variabilis, one may hypothesize that present Acidiferrobacteraceae couple sulfur oxidation to dissimilatory nitrite reduction, as recently described for Desulfurivibrio (44). However, S. variabilis, although isolated from a nitrate-reducing enrichment culture, was not able to grow using nitrate or nitrite as a terminal electron acceptor (42).

Given that S. variabilis also features a complete glycolysis pathway, tricarboxylic acid cycle, and respiratory chain, we speculate that Acidiferrobacteraceae at well H5-2 may participate in sulfate reduction (Fig. 6). However, given that we did not measure sulfate reduction and that genes involved in dissimilatory sulfate reduction/oxidation are commonly constitutive expressed (45, 46), the role of Acidiferrobacteraceae remains elusive.

Although autotrophic sulfate reduction is known, organic electron donors are preferentially used by sulfate-reducing bacteria (47). Recent work using stable and radiogenic carbon isotope analyses of PLFAs and carbon sources showed a preferred utilization of ancient carbon at H5-2 (48). The low Δ¹⁴C values for 10MeC_16:0, a lipid characteristic of some sulfate-reducing bacteria, suggests a strong involvement of present sulfate-reducing bacteria in the turnover of intermediates from organic matter breakdown.

We identified stable proportions of metagenome/metatranscriptome sequences linked to CAZymes (Fig. 5). The majority of these sequences were linked to cellulose, chitin, and pectin breakdown. A comparison of the overall proportions of CAZyme-affiliated sequences with values from other environments and ecosystems showed that the proportions determined for Hainich groundwater (1.0% to 1.7%) are at the lower end compared to values from temperate forest soils (1.6% to 2.1%) (49) but significantly lower in comparison to systems, such as paddy soil (3.3%% to 5.9%), that are accustomed to continuous exposure to complex polysaccharides (50).

In general, CAZyme-encoding genes constitute, on average, 3% of the genomes of most microorganisms (51). The uniformity of the encoded and expressed CAZyme potential did not match our expectations in seeing differences due to the diversity of carbon pools in terms of origins and ages along the Hainich groundwater transect. In contrast, the groundwater microbiome appears to be adjusted to a comparable stable carbon pool. Total organic carbon along the aquifer assemblages ranges from 2 to 4 mg liter⁻¹ (13). Easily utilizable carbon sources, eventually released from surface-derived input, are presumably quickly consumed along the flow path of the groundwater starting from the recharge area. As a result, the stable carbon pool is likely enriched with rather old carbon. Locally and temporally restricted spikes of easier utilizable carbon substrates can occur due to the recycling of necromass, which might be released from blooming populations that are eventually terminated by viruses according to the “kill-the-winner” model (52–54).

Despite their uniformity, encoded and expressed CAZyme functions were linked to different taxonomic groups dependent on the individual monitoring well (Fig. 5). It is apparent that the blooming of Bacteroidetes on the metatranscriptome level is a result of this group strongly contributing to complex polysaccharide breakdown, contrary to groups whose transcription levels were low although they possess pronounced degrading capabilities, such as Actinobacteria or Acidobacteria (51). As for overall taxonomic profiles (Fig. 2), Nitrospirae were partially highly abundant in CAZyme metagenome data sets (Fig. 5B), while their abundance in metatranscriptome data was reduced in comparison. Based on genomic profiling, the proportion of CAZyme-encoding genes in Nitrospirae is around 2% (51). Bacteroidetes are well known to feature comparable comprehensive enzymatic toolkits for complex polysaccharide breakdown (55). The number of CAZyme-encoding genes in Bacteroidetes genomes is often significantly higher than for other well-known polysaccharide-degrading organisms, for instance, Firmicutes (56). As a consequence, it is not surprising that Bacteroidetes mediate crucial roles in the context of polysaccharide breakdown in a multitude of different environmental settings, including but not limited to terrestrial environments, gut systems, as well as marine and freshwater systems (55, 57–59).

The use of metagenomics and metatranscriptomes in a complementary manner allowed us to holistically assess the importance of individual biogeochemical processes and taxonomic groups to microbiome functioning in the Hainich CZE. Being characterized by an overall nutrient scarcity, the Hainich CZE microbiome is governed by the spatially restricted availability of inorganic energy sources, resulting in hot spots for the chemolithoautotrophic processes nitrification and anaerobic ammonia oxidation. We could not confirm the importance of other biogeochemical processes linked to sulfur and iron cycling that were suggested to be relevant at individual monitoring wells based on hydrochemistry and PLFA profiles, due to the strong dominance of nitrogen cycling of encoded and expressed potentials. Nonetheless, based on our data, we hypothesize Acidiferrobacteraceae also to be involved in sulfate reduction, although this group was so far only associated with sulfur oxidation. The constitutive expression of polysaccharide-cleaving CAZymes highlighted the fact that microbial carbon assimilation by autotrophs is complemented by organic matter breakdown processes over the whole transect. Furthermore, Bacteroidetes appear to be key drivers for complex polysaccharide breakdown, also in oligotrophic groundwater.

The groundwater sampling sites are part of the Hainich CZE situated in the Hainich National Park in central Germany, a landscape characterized by sedimentary carbonate rocks and groundwater in karstified aquifers (12). A total of 11 groundwater wells accessing two superimposed aquifer assemblages (Hainich transect upper aquifer assemblage [HTU], Hainich transect lower aquifer assemblage [HTL]) were installed at five sites (H1 to H5) along a hillslope covering a distance of approximately 5.4 km, following the groundwater flow. Groundwater was reached at the following depths: H3-2, 22.5 m; H4-2, 12.7 m; H4-3, 12 m; H5-2, 65 m; and H5-3, 50 m at HTU; and H1-3, 5.1 m; H1-4, 7.5 m; H2-1, 31 m; H3-1, 47 m; H4-1, 58 m; and H5-1, 88 m at HTL. Oxygen availability is the major factor discriminating the two aquifer assemblages, with HTU being suboxic to anoxic, while HTL is oxic. The hydrochemistry of the two aquifer assemblages is strongly impacted by the limestone present and was thoroughly characterized previously (13). In short, the pH at all wells ranges from 7.1 to 7.2, while the redox potential varies between 200 and 400 mV. Multiparameter clustering suggested three biogeochemical zones in the transect, as follows: zone 1 (H1-3, H1-4, and H2-1 [cluster 1], H3-1, H3-2, and H4-1 [cluster 2], and H5-1 [cluster 3]), zone 2 (H4-2 and H4-3 [cluster 4]), and zone 3 (H5-2 and H5-3 [cluster 5]) (13). A more detailed description of the sampling sites can be found elsewhere (12, 13).

Sufficient biomass for subsequent nucleic extraction was collected by high-volume filtration, as described elsewhere (17), for sampling campaigns PNK66 and PNK69. The sampled wells included H3-2, H4-1, H4-2, H4-3, H5-1, and H5-2, and at each site, approximately 2,000 liters of groundwater was subjected to filtration through a glass fiber filter (0.3-μm pore size; Sterlitech, USA), with a flow rate of approximately 15 liters min^–1. Filters were immediately cooled at 4°C and transported to the laboratory, where they were stored at 80°C until being further processed for nucleic acid extraction. Replicated metagenome data sets and metatranscriptome data sets from pooled RNA preparations are available for PNK66. Complementary replicated metatranscriptome data sets were generated for PNK69 for wells H4-1, H5-1, H4-3, and H5-2.

Biomass-carrying filters were quartered, and one-quarter each was subjected to nucleic extraction. Filter pieces were washed in 15-ml conical tubes with 5 ml 2× AE buffer (10 mM Tris-HCl, 0.5 M EDTA [pH 9]) and 1.25 ml SDS (20% [vol/vol]) and subsequently vortexed for 1 min to facilitate detachment. After centrifugation (8,000 × g, 5 min), supernatants were transferred to new tubes and subsequently extracted with phenol-chloroform-isoamyl-alcohol (25:24:1) and chloroform-isoamyl-alcohol (24:1). The resulting supernatants were precipitated in fresh tubes by adding sodium acetate (3 M [pH 5.2]; final concentration, 0.3 M), glycogen (final concentration, 20 mg/ml), and two volumes of ethanol (100%), followed by 2 h of incubation at 4°C. The precipitated samples were spun down (8,000 × g, 5 min), and the obtained pellets were washed twice with ethanol (70%). Pellets were air-dried and resuspended in a suitable volume of AE buffer. Genomic DNA was quantified by spectrophotometry using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and its integrity checked by agarose gel electrophoresis. RNA was enriched from resuspended nucleic acids by digesting genomic DNA using Turbo DNase (Thermo Fisher, Schwerte, Germany), according to the manufacturer’s instructions.

Copy numbers of bacterial 16S rRNA genes were determined by quantitative PCR using an Mx3000P instrument (Agilent, Böblingen, Germany), Maxima SYBR green mastermix (Thermo Fisher Scientific, Germany), and the primer set Bac8Fmod/Bac338Rabc, as outlined before (60).

Library preparation and sequencing were done by LGC Genomics GmbH (Berlin, Germany). In brief, obtained genomic DNA or enriched mRNA was used, if possible, for technically replicated (three libraries per extract [per original filter]) sequencing library preparation. Aiming for a final insert size of 350 to 500 bp, sequencing libraries were generated using the Ovation rapid library system (NuGEN, San Carlos, CA, USA). Sequencing was carried out on an Illumina MiSeq platform in paired-end mode (2 × 300 bp). Per metagenome/metatranscriptome library, 1.7 to 3.0 million read pairs were obtained (Tables S3 and S4).

Sequence quality was assessed using FastQC (version 0.11.3; https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Adapter and quality trimming were performed with skewer (version 0.2.2 [61]) (parameters: –mode any –format sanger –mean-quality 20 –end-quality 20). Reads derived from rRNA genes were identified using SortMeRNA (version 2.1) (62) by database queries against stripped-down versions of SILVA (release 123) (63) and Rfam (release 12.2 [64]). Non-rRNA (gene)-derived metagenome/metatranscriptome reads were subsequently subjected to taxonomic and functional profiling. Detailed quality reports are available in the supplemental material (Supplementary Information 1 and 2 at https://osf.io/7432g/).

Metagenome/metatranscriptome reads were taxonomically assigned using Kaiju (65). Reads were translated into open reading frames, which were in turn used for string matching, aiming for maximum exact matches in reference sequences, applying a backward-search algorithm based on the one implemented in the Burrows-Wheeler transform (66, 67). NCBI nr (68) and bacterial and archaeal reference genomes deposited in NCBI RefSeq were used to construct a Kaiju protein database that was subsequently used for string matching. Kaiju was run in “greedy” mode, with a minimum score of 65 and allowing up to 5 amino acid substitutions. Taxonomic assignments were done based on retrieved database hits. In the case of multiple database hits of equal quality, a lowest common ancestor (LCA) approach was used to determine taxonomic affiliations. Given that the success rate of taxonomic annotation based on string matching is decreasing when looking at more resolved taxonomic levels (e.g., order or family), we complemented string matching with alignment-based taxonomy assignments. Reads were searched against NCBI RefSeq using DIAMOND (69) (applying an E value threshold of 0.001 and collecting up to 10 database hits) and annotated using the LCA algorithm implemented in MEGAN (version 6.0 [70]), with default settings. The applied settings for DIAMOND searches were the same as those outlined in the following section for functional annotation.

Metagenome/metatranscriptome reads were queried against a custom protein database constructed from available annotated prokaryotic genomes collected from NCBI RefSeq (68). Respective protein sequences were pooled and subsequently indexed for being used with DIAMOND (69). The DIAMOND searches were carried out as described above for taxonomic annotation. Search output was generated in the proprietary binary DAA format that facilitates data import into MEGAN (version 6.0 [70]). Functional annotations were done by parsing generated search output in MEGAN using available mapping files to the Kyoto Encyclopedia of Genomes and Genes (KEGG) (71). Sequence subsets of interest were eventually subsampled and subjected to in-depth analysis, for instance, determining taxonomic affiliations of functional subtranscriptomes.

CAZyme-affiliated sequences were determined by querying data sets against the dbCAN database (release 07202017) (72). Database queries were carried out with DIAMOND, applying the settings given above. Functional CAZyme modules of interest (cellulose, chitin, xylan, pectin turnover, and other hemicellulases) were defined by assembling related enzymatic functions based on enzyme commission numbers. A mapping file for functional annotation was created using all sequences and associated metadata (CAZy family, EC number) deposited in dbCAN. Functional annotation was facilitated by storing the mapping file as indexed sqlite database object and queried using custom python scripts (Supplementary Information 3 at https://osf.io/7432g/). The resulting annotations are based on matching dbCAN top hits for metagenome/metatranscriptome sequences queried against the formatted mapping file and defined CAZyme modules.

Statistical analyses were done with the R software framework (version 3.2.2 [73]) using the packages ape (version 3.3 [74]) and vegan (version 2.3.1 [75]), including their respective dependencies.

Figures were compiled using the R package ggplot2 (https://ggplot2.tidyverse.org/ [76]) and finalized with inkscape (https://inkscape.org/).

Sequence data were deposited at the European Nucleotide Archive (https://www.ebi.ac.uk/ena/) under BioProject number PRJEB28783.