Tree Species Shape Soil Bacterial Community Structure and Function in Temperate Deciduous Forests.

Amplicon-based analysis of 16S rRNA genes and transcripts was used to assess the effect of tree species composition on soil bacterial community structure and function in a temperate deciduous forest. Samples were collected from mono and mixed stands of Fagus sylvatica (beech), Carpinus betulus (hornbeam), Tilia sp. (lime), and Quercus sp. (oak) in spring, summer, and autumn. Soil bacterial community exhibited similar taxonomic composition at total (DNA-based) and potentially active community (RNA-based) level, with fewer taxa present at active community level. Members of Rhizobiales dominated at both total and active bacterial community level, followed by members of Acidobacteriales, Solibacterales, Rhodospirillales, and Xanthomonadales. Bacterial communities at total and active community level showed a significant positive correlation with tree species identity (mono stands) and to a lesser extent with tree species richness (mixed stands). Approximately 58 and 64% of indicator operational taxonomic units (OTUs) showed significant association with only one mono stand at total and active community level, respectively, indicating a strong impact of tree species on soil bacterial community composition. Soil C/N ratio, pH, and P content similarly exhibited a significant positive correlation with soil bacterial communities, which was attributed to direct and indirect effects of forest stands. Seasonality was the strongest driver of predicted metabolic functions related to C fixation and degradation, and N metabolism. Carbon and nitrogen metabolic processes were significantly abundant in spring, while C degradation gene abundances increased from summer to autumn, corresponding to increased litterfall and decomposition. The results revealed that in a spatially homogenous forest soil, tree species diversity and richness are dominant drivers of structure and composition in soil bacterial communities.

Introduction

Forests offer many ecosystem services of ecological and economic significance, provide a diverse, multi-layered habitat for most terrestrial plants and animals, and resources for pan class="Species">humans ranging from timber to recreational facilities (Cardenas et al., 2015; Wood et al., 2017). Forests are central hubs of class="Chemical">primary class="Chemical">productivity and a major sink of class="Chemical">pan class="Chemical">carbon, due to the large concentration of trees and other vegetation (Bonan, 2008). Thus, the resulting interplay of ecosystem functions within forests is dynamic and often performed by both biotic and abiotic forest components.

In particular, microorganisclass="Disease">ms facilitate the exchange of comclass="Chemical">pounds between forest vegetation and forest soils. Forest soils are essential biological matrices in which microbial communities execute key ecosystem functions, such as biogeochemical cycling, through decomclass="Chemical">position and mineralization class="Chemical">processes mediated by class="Chemical">prokaryotes and fungi (Uroz et al., 2016a; Llclass="Chemical">pan class="Disease">adó et al., 2017). Factors influencing soil microbial community structure are crucial for predicting how bacteria-mediated processes drive ecosystem responses to environmental changes (Nemergut et al., 2014). Soil bacterial communities are shaped by several edaphic factors, including soil texture and chemistry, as well as biotic factors such as plant roots and associated mycorrhizal activity, aboveground litter, and other decomposing organic matter (Thoms et al., 2010; Lladó et al., 2018).

Compared to grassland or agricultural ecosysteclass="Disease">ms, only a limited number of studies are available on the effects of forest tree sclass="Chemical">pecies on bacterial community structure and function. These studies revealed comclass="Chemical">plex bacterial community organization that results from dynamic forest resources and class="Chemical">processes. Driven by differences in soil class="Chemical">pH, beech stands exhibited a higher bacterial diversity comclass="Chemical">pared to sclass="Chemical">pruce stands (Nacke et al., 2011). In mesocosm exclass="Chemical">periments comclass="Chemical">prising five mono and mixed sclass="Chemical">pecies stands, soil bacterial community richness and evenness were strongly influenced by the tree sclass="Chemical">pecies beech and ash (Pfeiffer et al., 2013). A class="Chemical">pan-class="Chemical">pan class="Species">European study on land-use intensity and microbial co-occurrence indicated that forest soil communities form far richer networks than those of grassland or farm soils (Creamer et al., 2016). Although soil characteristics (particularly pH) have been frequently reported as strong drivers of microbial diversity (Lauber et al., 2008; Kaiser et al., 2016), tree species have been shown to exhibit a stronger impact on community structure than the soil environment (Bonito et al., 2014).

Forest trees exert a stronger influence on soil than other perennial vegetation due to their longevity and the lack of soil management in forests. Over time, changing features, such as forest canopy, root biomass, root exudates, and class="Chemical">oxygen and class="Chemical">pan class="Chemical">water consumption, change soil temperature, chemistry, porosity, and soil moisture (Augusto et al., 2002, 2015), which in turn affect belowground bacterial communities (Uroz et al., 2016b). Seasonal climate affects tree species-dependent patterns of leaf phenology, leaf litter, and dead wood quality, as well as root biomass development (Goodale et al., 2015; Lladó et al., 2017; Žifčáková et al., 2017). Consequently, belowground carbon (C) allocation, available nitrogen (N), carbon/nitrogen ratio (C/N), and other nutrients of belowground bacterial communities also vary seasonally. This results in selective enrichment of certain taxonomic groups and metabolic processes (Rasche et al., 2011; Yokobe et al., 2018). Tree species such as beech (Fagus), spruce (Picea), and pine (Pinus) have demonstrated stronger acidifying effects on soils due to slow decomposition rates than ash, hornbeam, and lime (Augusto et al., 2002; Uroz et al., 2016b). Soil-acidifying tree species promoted the presence of a high proportion of acidophilic bacterial taxa (Uroz et al., 2016b). Members of Rhizobiales and Burkholderiales orders were influenced by pH differences arising from root exudates of Acer, Betula, Fagus, and Quercus spp. rather than by soil properties or geographic distance (Landesman et al., 2014). The degree of influence on soil bacterial community structure and diversity depends on the tree species, stand type, and spatial arrangement, i.e., mono species or mixed species stands (Klimek et al., 2016; Uroz et al., 2016a). Plant root exudates, including those from trees, produce a highly variable quality and quantity of carbohydrates, organic acids, and other signaling molecules depending on different factors, including tree species, age, as well as environmental and edaphic conditions (Vranova et al., 2013; Urbanová et al., 2015; Zhalnina et al., 2018). Subsequently, root exudates were shown to promote the abundance of taxa beneficial to specific plant species (Zhalnina et al., 2018). Thus, forest trees affect soil bacterial composition and function through direct and indirect processes.

To our knowledge, field studies reporting how different tree species and their spatial arrangement affect soil bacterial communities in deciduous forests are not available. The few studies that focused on how forest stands influence soil bacterial community structure at both entire and active bacterial community level, were conducted in soils with substantial spatial heterogeneity (Chodak et al., 2016; Kpan class="Species">limek et al., 2016; Siles and Margesin, 2017; Curd et al., 2018). Although forest trees demonstrated sclass="Chemical">pecies-declass="Chemical">pendent selection of sclass="Chemical">pecific bacterial taxa, the imclass="Chemical">pact of mono and mixed stands in sclass="Chemical">patially homogeneous soil, i.e., soil showing uniform class="Chemical">physical comclass="Chemical">position and constant elevation, on soil bacterial community structure was not exclass="Chemical">plored (Uroz et al., 2016b).

The aim of this study was to investigate the influence of tree species on soil bacterial communities in an old-growth bropan class="Disease">ad-leaved temclass="Chemical">perate forest, the Hainich National Park. We hyclass="Chemical">pothesized that (1) tree sclass="Chemical">pecies identity drives bacterial community diversity and structure, resulting in taxonomic differences at entire and class="Chemical">potentially active bacterial community level, and (2) metabolic functions are also driven by tree sclass="Chemical">pecies, and follow stand-sclass="Chemical">pecific comclass="Chemical">position. We class="Chemical">predicted that (3) class="Chemical">pure (mono sclass="Chemical">pecies) and mixed stands are resclass="Chemical">ponsible for changes in soil class="Chemical">physicochemical class="Chemical">proclass="Chemical">perties, and that bacterial communities in different stands are shaclass="Chemical">ped by season but to a lesser extent than by tree sclass="Chemical">pecies. We evaluated the hyclass="Chemical">potheses by comclass="Chemical">paring 16S rRNA gene and transcriclass="Chemical">pt amclass="Chemical">plicon-derived soil bacterial community data, and functional class="Chemical">predictions with soil class="Chemical">physicochemical class="Chemical">parameters of soils from mono and mixed deciduous tree stands.

Materials and Methods

Sample Site Description

The study was conducted over a 25-ha area in the Hainich National Park located in Thuringia, Germany (Figures 1A,B). The park is the largest unmanaged deciduous, broclass="Disease">ad-leaved forest ecosystem in central Germany, situated near the village of Weberstedt (350 m a.s.l; 51°05′37.0 N, 10°30′10.6 E) (Mölder et al., 2006; Schmidt et al., 2015). The samclass="Chemical">pling area lies over uniform bedrock of Triassic class="Chemical">pan class="Species">limestone covered by up to 50 cm of Loess and soil with silt loam/silty clay loam Cambisol (Schmidt et al., 2015). Since the 1960s, the area underwent minor management activity, which was discontinued after its declaration as a national park in 1997 and is therefore described as a natural ecosystem (Kubisch et al., 2015). Core tree species in the mixed-forest ecosystem include the European beech (Fagus sylvatica), small-leaved and large-leaved lime (Tilia cordata Mill. and T. platyphyllos Scop., respectively), oak (Quercus petraea and Q. robur), and hornbeam (Carpinus betulus) (Mölder et al., 2006). Sample sites were selected based on the following criteria: (1) each stand should have a similar surrounding tree species composition displayed, (2) trees must be actively growing and have a well-defined canopy, and (3) homogeneity in soil parameters including color, texture, drainage, and slope of sampling area at initial sampling (Schmidt et al., 2015). Individual stands consisted of four to eight trees, each with an area between 68 and 313 m2. Mono-species stands (hereafter referred to as mono stands) comprised beech, hornbeam, lime, and oak, while mixed-species stands (hereafter referred to as mixed stands) consisted of three out of the four tree species: BHL (beech, hornbeam, lime), BHO (beech, hornbeam, oak), BOL (beech, oak, lime), and HOL (hornbeam, oak, lime). Each stand was replicated six times, resulting in a total of 48 plots (Figure 1C).

Figure 1

Area of study in the Hainich national park. (A) Boundary of the Hainich national park (51°05′37.0 N, 10°30′10.6 E). (B) Map of Germany showing location of the Hainich National Park in Thuringia, Germany. (C) Sampling area showing distribution of mono stands and mixed stands used for sample collection. Abbreviations stand for the following tree species associations in mixed stands: beech-hornbeam-lime (BHL), beech-oak-hornbeam (BHO), beech-oak-lime (BOL), and hornbeam-oak-lime (HOL). The map was generated by using data from the following data providers: Esri, HERE, DeLorme, Intermap, increment P. Corp., GEBCO, USGS, FAO, NPS, NRCAN, GeoBase, IGN, Kadaster, NL, Ordnance Survey, Esri Japan, METI, Esri China (Hong Kong), swisstopo, Mapmyindia, @OpenStreatMap contributors, and the GIS user community.

Area of study in the Hainich national park. (A) Boundary of the Hainich national park (51°05′37.0 N, 10°30′10.6 E). (B) Map of Germany showing location of the Hainich National Park in Thuringia, Germany. (C) Sampling area showing distribution of mono stands and mixed stands used for sample collection. Abbreviations stand for the following tree species associations in mixed stands: beech-hornbeam-class="Species">lime (BHL), beech-oak-hornbeam (BHO), beech-oak-class="Chemical">pan class="Species">lime (BOL), and hornbeam-oak-lime (HOL). The map was generated by using data from the following data providers: Esri, HERE, DeLorme, Intermap, increment P. Corp., GEBCO, USGS, FAO, NPS, NRCAN, GeoBase, IGN, Kadaster, NL, Ordnance Survey, Esri Japan, METI, Esri China (Hong Kong), swisstopo, Mapmyindia, @OpenStreatMap contributors, and the GIS user community.

Sampling and Environmental Nucleic Acid Isolation

Soil samples were collected in spring (April), summer (July), and autumn (September) of 2012. Two soil cores (10 cm in diameter, 5 cm depth) from the A horizon (topsoil) were extracted from bulk soil at randomly selected polar sites within a plot and pooled to generate a composite sample. Determinations of soil physicochemical properties such as class="Chemical">nitrogen (N), class="Chemical">pan class="Chemical">carbon (C), phosphorus (P) content, carbon/nitrogen (C/N) ratio and moisture content were performed and described by Schmidt et al. (2015) and are presented in Supplementary Figure S1.

Extraction of envclass="Chemical">ironmental DNA and RNA was class="Chemical">performed with the MoBio Powersoil DNA isolation kit and RNA Powersoil Total RNA isolation kit, resclass="Chemical">pectively, as recommended by the manufacturer (MoBio Laboratories, Carlsbclass="Chemical">pan class="Disease">ad, USA). Minor modifications were performed in the extraction of total RNA as follows: after eluting bound RNA from the RNA capture column, a 2.5-fold volume of ethanol and 0.1-fold volume of 3 M sodium acetate (pH 5.2) were added to the RNA, vortexed, and stored at −80°C until required. Solutions were centrifuged at 10,000 × g for 1 h (4°C) to recover RNA. Subsequent ethanol/sodium acetate supernatants were discarded, and RNA pellets were dried before suspension in 100 μl distilled water. Residual DNA was removed from extracted RNA by using TURBO DNA-Free kit (Ambion Applied Biosystems, Darmstadt, Germany). Absence of DNA was confirmed according to Wemheuer et al. (2012). Resulting RNA was concentrated with the Rneasy MinElute cleanup kit (Qiagen GmbH, Hilden, Germany). Quantification of DNA and RNA was performed with the NanoDrop ND-1000 UV-Vis spectrophotometer, following the instructions of the manufacturer (Peqlab Biotechnologie GmbH, Erlangen, Germany). Purified RNA (approximately 300–400 ng) was converted to cDNA using SuperScript III reverse transcriptase (Invitrogen, Karlsruhe, Germany) and the bacterial reverse primer V5rev_B 5′-CTATGCGCCTTGCCAGCCCGCTCAG-MID-CCGTCAATTCMTTTGAGT-3′ (Wang and Qian, 2009).

Envclass="Chemical">ironmental DNA and cDNA were used as temclass="Chemical">plates to amclass="Chemical">plify the V3–V5 regions of the 16S rRNA gene by PCR. The 50 μl PCR reaction mixture contained 25 ng of envclass="Chemical">pan class="Chemical">ironmental DNA or cDNA as template, 10 μl of 5X Phusion GC buffer, 0.2 μM of each of the four deoxynucleoside triphosphates (dNTPs), 0.4 μM of each primer, 2.5 μl DMSO, 0.15 μl of 25 mM MgCl2, and 1 U of Phusion high-fidelity DNA polymerase. The V3–V5 region was amplified with the following set of primers comprising the Roche 454 pyrosequencing adaptors (underlined), a key (TCAG), a unique 10-bp multiplex identifier (MID), and template-specific sequence per sample: the forward primer V3for_B (5′-CGTATCGCCTCCCTCGCGCCATCAG-MID-TACGGRAGGCAGCAG-3′), and reverse primer V5rev_B (see above) from Liu et al. (2007) and Wang and Qian (2009), respectively. Positive and negative control reactions were prepared by using bacterial DNA or substituting water (DNA and RNA-free) instead of the environmental nucleic acids. All PCR reactions were performed in triplicate employing the following thermal cycling parameters for amplification: initial denaturation at 98°C for 5 min, followed by 25 cycles of denaturation at 98°C for 45 s, annealing at 65°C for 45 s and extension at 72°C for 30 s, and a final extension at 72°C for 5 min. Resulting amplicons were analyzed by gel electrophoresis, pooled, and purified with the Qiagen Qiaquick gel extraction kit as recommended by the manufacturer (Qiagen GmbH). Quantification of amplicons was determined by using the Quant-iT dsDNA BR assay kit and Qubit fluorometer as recommended by the manufacturer (Invitrogen GmbH, Karlsruhe, Germany). Amplicon sequencing was performed by the Göttingen Genomics Laboratory with the 454 GS-FLX+ pyrosequencer and titanium chemistry as recommended by the manufacturer (Roche, Mannheim, Germany).

Processing of 16S rRNA Gene Sequence Data and Statistical Analyses

The resulting 16S rRNA gene and transcript sequences were processed and analyzed with the QIIME (1.9.1) software package (Caporaso et al., 2010) by employing the scheme outlined by Schneider et al. (2015). This involved removal of sequences shorter than 300 bp, containing unresolved nucleotides, exhibiting a low average quality score (<25) or harboring long homopolymers (>8 bp). Forward and reverse primer sequences were removed with the split_libraries.py script. Remaining reverse primer sequences, pyrosequencing noise and chimeric sequences, were removed with cutclass="Disease">adaclass="Chemical">pt (Martin, 2011), Acacia (Bragg et al., 2012), and UCHIME (Edgar et al., 2011), resclass="Chemical">pectively. Oclass="Chemical">perational taxonomic units (OTUs) were assigned at 97% genetic similarity (sclass="Chemical">pecies-level) by emclass="Chemical">ploying the class="Chemical">pick_oclass="Chemical">pen_reference_otus.class="Chemical">py scriclass="Chemical">pt against the class="Chemical">pan class="Disease">Silva database (release 128) (Quast et al., 2013). Taxonomic classification of OTUs was performed by parallel_assign_taxonomy_blast.py script the Silva SSU database release 128. The filter_otu_table.py script was used to remove singletons, chloroplast sequences, extrinsic domain OTUs, and unclassified OTUs.

For all statistical tests, a p ≤ 0.05 was regarded as significant. Alpha and beta diversity indices and rarefaction curves were calculated with QIIME by using alpha_rarefaction.py employing the same level of surveying effort (13,000 randomly selected bacterial sequences per sample). class="Disease">Additional analyses were class="Chemical">performed in R (R Develoclass="Chemical">pment Core Team, 2017). Non-metric multidimensional scaling (NMDS) was class="Chemical">performed with the “vegan” class="Chemical">package (Oksanen et al., 2019), emclass="Chemical">ploying weighted UniFrac distance matrix calculated by beta_diversity_through_class="Chemical">plots.class="Chemical">py (Lozuclass="Chemical">pone et al., 2011). Envclass="Chemical">pan class="Chemical">ironmental parameters were fitted on the NMDS with the envfit function of the “vegan” package in R (Gergs and Rothhaupt, 2015). A combination of analysis of similarities (ANOSIM) and permutational multivariate analysis of variance (PERMANOVA), implemented in QIIME through compare_categories.py were used to improve robustness of multivariate analyses on the effect of tree stands on bacterial community using weighted UniFrac distance matrix (Hartmann et al., 2015). Mantel test was performed using the “vegan” package in R to test correlations between soil physicochemical parameters and soil bacterial communities. Association networks between tree species and OTUs were determined by mapping significant point biserial correlation values calculated by the “indicspecies” package in R (De Cáceres, 2013). Subsequent network visualizations for taxa/tree stand associations were generated with Cytoscape v3.5 by using the “edge-weighted spring embedded layout” algorithm, whereby network edges were weighted by association value (Shannon et al., 2003; Cline et al., 2007). Prediction of functional pathways and corresponding enzymes was performed with Tax4Fun (Aßhauer et al., 2015). Subsequent NMDS and boxplots were calculated with the vegdist function in “vegan” and “ggplot2” package, respectively (Gergs and Rothhaupt, 2015; Wickham, 2016).

Accession Numbers

The 16S rRNA gene and transcript sequences were deposited in the National Center for Biotechnology Information (NCBI) Sequence Repan class="Disease">ad Archive (SRA) under the accession number PRJNA342582.

Results and Discussion

Edaphic Properties Reflect Stand-Related Variations

Soil properties displayed significantly different stand-specific variations (p < 0.05) except for moisture content (Supplementary Figure S1). class="Chemical">Carbon content only differed in hornbeam and the HOL mixed stand. As class="Chemical">pan class="Species">European forests store large stocks of soil organic carbon (De Vos et al., 2015), soil carbon was not expected to show significant stand-specific concentrations. Among mono stands, lime exhibited the highest mean pH (5.9 ± 0.6). Beech and oak shared a low pH environment (4.6 ± 0.3 and 4.5 ± 0.5, respectively). Mixed stands that included lime displayed higher pH values (BOL 5.4 ± 0.4 and HOL 5.2 ± 0.5) than mixed stands with beech and hornbeam (BHL 5 ± 0.4 and BHO 4.8 ± 0.8).

C/N ratios between 8 and 16 indicate complete breakdown of organic material and, consequently, higher class="Chemical">nitrogen availability while class="Chemical">pan class="Chemical">phosphorus (P) content, as an essential component for plant growth, indicates organic matter richness and quality (Lauber et al., 2008; Lang et al., 2016). Therefore, we used C/N and P content as indicators of soil fertility (Supplementary Figure S1). Beech and oak mono stands showed high C/N ratios (15.2 ± 0.7 and 16.1 ± 1.3, respectively) compared to lime mono stands (12.4 ± 0.6). Correspondingly, beech and oak mixed stands (BHO and BOL) also exhibited higher C/N ratios compared to lime and hornbeam (BHL) mixed stands (Supplementary Figure S1). Plant available P was consistently higher in soils of hornbeam (585 ± 240 mg/kg) and lime mono stands (536 ± 109 mg/kg) compared to beech (340 ± 41 mg/kg) and oak (335 ± 58 mg/kg). Furthermore, mixed stands with hornbeam and lime (BHL and HOL) exhibited significantly higher P content than mixed stands with beech and oak (BHO and BOL).

Low C/N ratios, high available N and P, as well as high pH, as observed in class="Species">lime mono stands, class="Chemical">promoted tree class="Chemical">productivity (Schmidt et al., 2015). In contrast, beech and oak mono stands exhibited higher C/N ratios, low soil available class="Chemical">pan class="Chemical">nitrogen and low pH. An explanation could be the decomposition rates of leaf litter. In both mono and mixed stands, lime and hornbeam leaf litter decomposed faster compared to beech and oak leaf litter (Schmidt et al., 2016). Consequently, faster decomposition rates contributed to the observed effect on soil chemistry by a faster release of nutrients (Jacob et al., 2009; Scheibe et al., 2015; Schmidt et al., 2016).

Soil characteristics in mono stands were also recorded in associated mixed stands. This indicates that tree identity and stand type may give rise to distinct microenvclass="Chemical">ironments, whereby mixed stands resulted in intermediate effects in stand class="Chemical">productivity comclass="Chemical">pared to mono stands (Schmidt et al., 2015). The results showed a stronger influence of tree sclass="Chemical">pecies identity over tree diversity on soil characteristics. Our results are suclass="Chemical">pclass="Chemical">ported by similar findings within the Hainich region, which consistently showed a more acidic but class="Chemical">potentially nutrient limited soil envclass="Chemical">pan class="Chemical">ironment in beech- or oak-dominated stands, compared to lime- or hornbeam-dominated stands (Brunet et al., 1997; Falkengren-Grerup et al., 1998; Salehi et al., 2007; Langenbruch et al., 2012; Berger and Berger, 2014).

Stand-Specific Patterns in Bacterial Richness and Diversity

Soil bacterial community structure was determined based on 16S rRNA gene and transcript analyses. After processing, quality-filtering and taxonomic clustering at 97% similarity, we recovered 40,385 for entire (DNA-based) and 52,277 OTUs for active (RNA-based) bacterial communities (Supplementary Table S1). The higher number of bacterial OTUs in the active community is in contrast to some other studies of bacterial community diversity in forest (Baldrian et al., 2012; Romanowicz et al., 2016) and grassland soils (Herzog et al., 2015). However, Gill et al. (2017) reported higher RNA-based OTU counts from urban soils, similar to our results. We attribute the different abundances in total and active communities to the detection of very low-abundant but active rare taxa at RNA level.

Diversity indices showed several consistent patterns in mono and mixed tree stands with respect to species richness and evenness, as indicated by Shannon (H’), Chao 1, and phylogenetic diversity (PD) (p < 0.05; Supplementary Figure S2). Rarefaction curves were close to saturation, and Good’s coverage estimator across all stands remained above 75%, indicating that the sampling effort recovered most of the bacterial diversity.

Among mono stands, class="Species">lime mono stands showed the highest soil bacterial richness (H′ 10.5 ± 0.3; PD 169 ± 33) at total and active community level. Similarly, among mixed stands (H′ 10.4 ± 0.3; PD 170 ± 32) soil bacterial richness was the highest in class="Chemical">pan class="Species">lime-associated mixed stands (BOL and HOL). Oak mono stands showed significantly lower bacterial richness and diversity (H′ 9.8 ± 0.4; PD 140 ± 27) compared to other mono stands. We attributed differences in bacterial richness and diversity to tree species-specific effects on the bacterial community. Changes in bulk soil chemistry, i.e., pH differences and nutrient availability, drive bacterial community structure and diversity (Jeanbille et al., 2016; Kaiser et al., 2016). These changes are a direct result of leaf litter and root exudates in soil, which subsequently impact belowground soil bacterial communities (Thoms et al., 2010). In our study, lime mono stands exhibited the highest bacterial richness and diversity and oak the lowest. The favorable physiological conditions (higher pH, low C/N ratio, high exchangeable P) in lime and hornbeam mono stands promoted higher bacterial diversity, compared to the acidic environment of beech and oak mono stands. This is supported by studies in a mixed deciduous forest, which showed similar trends in bacterial community diversity in mono stands and two-species mixed stands of lime, oak and alder (Alnus sp.), birch (Betula sp.), larch (Larix sp.), and spruce (Picea sp.) (Šnajdr et al., 2013; Urbanová et al., 2015). Supplementary Figure S3 provides further details on alpha diversity metrics of OTUs at total community (DNA-derived) and potentially active (RNA-derived) levels.

Tree Species Effect on Bacterial Community Structure

Multivariate analyses confirmed that soil bacterial communities can be differentiated according to tree species (tree species, p = 0.001, r 2 = 0.506) (Figure 2 and Supplementary Table S2). Bacterial communities of each tree species clustered together at total and active community level (template, p = 0.001, r 2 = 0.439). We furthermore detected a pattern in which communities in mono stands separated from those in mixed stands (stand type, p = 0.043, r 2 = 0.068). Exceptions were bacterial communities in hornbeam mono stands, which grouped with mixed stands. Results of the NMDS were supported by ANOSIM and PERMANOVA tests, which indicated that tree species was a strong driver of bacterial community structure at total and active community level (p = 0.001, ANOSIM and PERMANOVA; Supplementary Table S2). Stand type also exhibited a strong influence on community structure at the potentially active community level (p = 0.002 and p = 0.044, ANOSIM and PERMANOVA, respectively). The results supported the trends observed in alpha diversity analyses, in which bacterial richness and evenness show stand-specific variation. Thus, these results confirm hypothesis (1), that tree species drive bacterial community composition.

Figure 2

Non-metric multidimensional scaling (NMDS) analysis of soil bacterial communities in mono and mixed stands. (A) NMDS plot for total and active community. (B,C) NMDS plots for total and active community, respectively, showing the goodness of fit for soil environmental parameters. Ordination employed weighted UniFrac distance matrices of OTUs. Samples are grouped according to season, whereby each data point represents the summarized taxonomic data of six replicate stands. Arrows outside the plot area have no statistical significance and are only a visual aid of observed trends. C/N ratio (C/N), phosphorus content (P), nitrogen content (N), soil moisture (MC), operational taxonomic unit (OTU).

Non-metric multidimensional scaling (NMDS) analysis of soil bacterial communities in mono and mixed stands. (A) NMDS plot for total and active community. (B,C) NMDS plots for total and active community, respectively, showing the goodness of fit for soil envclass="Chemical">ironmental class="Chemical">parameters. Ordination emclass="Chemical">ployed weighted UniFrac distance matrices of OTUs. Samclass="Chemical">ples are grouclass="Chemical">ped according to season, whereby each data class="Chemical">point reclass="Chemical">presents the summarized taxonomic data of six reclass="Chemical">plicate stands. Arrows outside the class="Chemical">plot area have no statistical significance and are only a visual aid of observed trends. C/N ratio (C/N), class="Chemical">pan class="Chemical">phosphorus content (P), nitrogen content (N), soil moisture (MC), operational taxonomic unit (OTU).

Seasonal effects in our data were less pronounced (season, p = 0.858, r 2 = 0.0142) (Figure 2A, Supplementary Table S2), although a loose clustering of communities was observed (Figures 2B,C). ANOSIM and PERMANOVA analyses indicated that season was significantly correlated with bacterial communities at total community level (p = 0.464 and 0.529, respectively), but ANOSIM showed a significant seasonal correlation at active community level (p = 0.038). Although temporal seasonality is an important factor in temperate deciduous forests, which influences plant phenology, and subsequent root exudation and nutrient uptake processes (Oh et al., 2012), we did not observe a strong link between season and bacterial taxonomic composition in this study. Correspondingly, a warming experiment of a temperate forest soil only showed a shift in bacterial community composition after 20 years and only in the organic horizon (DeAngelis et al., 2015). Seasonal effects on soil bacteria at active community level can be explained by small changes driven by “conditionally rare taxa” (CRT), which periodically increase envclass="Chemical">ironmental during favorable conditions (Shclass="Chemical">pan class="Disease">ade et al., 2014). Moreover, CRT may maintain low abundance within the environment but during fluctuating conditions, as during seasonal changes, react by rapidly increasing the abundance of the ribosomal content (Hausmann et al., 2016). The active bacterial community is more sensitive to environmental disturbance than the total community (Pochon et al., 2017), as environmental RNA almost exclusively originates from viable organisms (Moran et al., 2013; Pochon et al., 2017). These effects are likely to produce distinct seasonal communities at active community level.

Edaphic parameters also contributed significantly in shaping the bacterial community (p = 0.001; Figures 2B,C; Supplementary Table S2). The Mantel statistic showed a significant positive correlation of soil parameters with bacterial communities (Supplementary Table S2). In particular, a strong influence of soil pH, P, N, and soil moisture was recorded in class="Species">lime and hornbeam mono stands, and in the corresclass="Chemical">ponding BHL, BOL, and HOL mixed stands comclass="Chemical">pared to the other stands. Bacterial communities in beech and oak mono stands and the BHO mixed stand resclass="Chemical">ponded to an increase in C/N ratio. Nutrient availability in soil has been linked to soil bacterial community structure (Bergkemclass="Chemical">per et al., 2016). Our results show that the nutrient-rich soils of class="Chemical">pan class="Species">lime and hornbeam mono stands possess bacterial community structures, which are distinct from those in comparatively nutrient-reduced soils in beech and oak mono stands.

Rather than tree species richness (mixed stands), tree species identity (mono stands) contributes to bacterial diversity and structure, as tree species showed a stronger positive correlation with bacterial communities than stand type (tree species r 2 = 0.506; stand type r 2 = 0.068). A study in a mixed deciduous forest demonstrated a similar correlation of microbial diversity with tree species identity, while species richness influenced microbial composition (Khlifa et al., 2017). Plant-specific traits such as fine root biomass and density, and leaf litter quality contributed to nutrient availability in soil (Thoclass="Disease">ms et al., 2010). Beech stands were shown to have higher fine root biomass than hornbeam or class="Chemical">pan class="Species">lime (Jacob et al., 2013), which could have a stronger impact on bacterial communities than other mono stands. However, beech showed no overyielding of root biomass in mixed stands (Langenbruch et al., 2012; Jacob et al., 2013). This explains that the impact of mixed stands on bacterial communities is often similar to the constituent mono stands, as little competition exists in fine roots. Furthermore, a significant change in soil bacterial community was observed with increasing horizontal distance from trees and demonstrated shifts in nutrient availability, such as decreasing concentrations of organic nitrogen compounds (Nacke et al., 2016). These changes correspond to a decrease in species-dependent root biomass density and subsequent nutrient availability.

Root exudation is plant-species dependent (Zhalnina et al., 2018). The release of exudates into the soil can be direct or regulated further through interactions of arbuscular mycorrhizal or ectomycorrhizal fungi with plant roots (Meier et al., 2013; Kaiser et al., 2015; Canarini et al., 2019). Current understanding of the transfer of root exudates from the rhizosphere to the bulk soil is still poorly understood; however, evidence suggests that nutrient transfer by ectomycorrhiza into the soil results in rapid changes in nutrient conditions of the surrounding soil (Gorka et al., 2019). All tree species in our study show ectomycorrhizal colonization (Liese et al., 2018), and the degree of mycorrhiza influence on bacterial community structure is also likely dependent on tree species.

General Patterns in Bacterial Community Composition and Structure

The entire dataset comprised 40 bacterial phyla, 155 classes, 385 orders, 704 families, and 1,552 genera. Dominant phyla (relative abundance > 1%) across the entire dataset comprised Acidobacteria, Proteobacteria, Actinobacteria, Bacteroidetes, Gemmatimonpan class="Disease">adetes, Chloroflexi, Latescimicrobia, and Nitrosclass="Chemical">pirae (Figure 3). Acidobacteria dominated at total bacterial community level (37.2% at DNA level and 23.2% at RNA level) and Proteobacteria at active community level (36.8% at DNA level and 56% at RNA level). Generally, soils contain a few dominant bacterial class="Chemical">phyla, mainly Acidobacteria, Actinobacteria, Bacteroidetes, and Proteobacteria (Fierer et al., 2012) that have been found consistently in both cultivated and natural envclass="Chemical">pan class="Chemical">ironments (Shi et al., 2011; Shange et al., 2012; Ferrenberg et al., 2013; Shen et al., 2014; Schneider et al., 2015).

Figure 3

Mean relative abundance of soil bacterial communities in mono and mixed stands. Results show bacterial community composition at order level, whereby each stand represents the average abundance of n = 3 seasons (spring, summer, and autumn). Total and potentially active taxa were derived from environmental DNA and RNA, respectively. Taxa with a relative abundance less than 1% across all stands were grouped into “rare bacterial orders.” Abbreviations for mixed stands: beech-hornbeam-lime (BHL), beech-hornbeam-oak (BHO), beech-oak-lime (BOL), hornbeam-oak-lime (HOL).

Mean relative abundance of soil bacterial communities in mono and mixed stands. Results show bacterial community composition at order level, whereby each stand represents the average abundance of n = 3 seasons (spring, summer, and autumn). Total and potentially active taxa were derived from envclass="Chemical">ironmental DNA and RNA, resclass="Chemical">pectively. Taxa with a relative abundance less than 1% across all stands were grouclass="Chemical">ped into “rare bacterial orders.” Abbreviations for mixed stands: beech-hornbeam-class="Chemical">pan class="Species">lime (BHL), beech-hornbeam-oak (BHO), beech-oak-lime (BOL), hornbeam-oak-lime (HOL).

Bacterial community composition was stand-specific (Supplementary Figures S4, S5) and showed dependence on pH and nutrient availability (Supplementary Figure S1). Copiotrophic members, such as Alpha- and Betaproteobacteria, Actinobacteria, and Bacteroidetes, are class="Disease">adaclass="Chemical">pted to nutrient-rich envclass="Chemical">pan class="Chemical">ironments, as found in the soils of lime stands. In contrast, members of Acidobacteria are recognized as oligotrophs, which are adapted to nutrient-limiting conditions (Fierer et al., 2012; Shange et al., 2012; Koyama et al., 2014). Beech stands, which have been frequently associated with nutrient-limited acidic soils, favored acidobacterial groups (Uroz et al., 2011; Jeanbille et al., 2016; Lladó et al., 2016; Kielak et al., 2016b; Colin et al., 2017).

Proteobacteria were evenly distributed across all sites with slight differences between mono stands and mixed stands, at total community level (Figure 3 and Supplementary Table S3). The Rhizobiales order within Alphaproteobacteria were more abundant in class="Species">lime and hornbeam mono stands at total community and active community level than in beech and oak mono stands (Figure 3 and Suclass="Chemical">pclass="Chemical">plementary Table S3). A similar trend was observed in the mixed stands BHL and HOL comclass="Chemical">pared to BHO and BOL. At genus level, we observed that these differences were due to high abundances of unidentified members of the Xanthobacteraceae family, Brclass="Chemical">pan class="Disease">adyrhizobium, Rhizomicrobium, and Variibacter. Rhizobiales are broadly associated with nitrogen fixation, plant pathogenicity, and organic matter decomposition (Carvalho et al., 2010). The Rhodospirillales order, which comprised mostly uncultured groups at genus level, showed high abundance in soils of beech and lime mono stands, due to the presence of the Acetobacteraceae family, in beech mono stands (up to 4.7% at DNA and RNA level) and uncultured Rhodospirillum DA 111 in beech and lime mono stands (up to 11.3% at DNA and RNA level). Acetobacteraceae contains some acidophilic genera of acetic acid bacteria (AAB), which are adapted to acidic environments (Mamlouk and Gullo, 2013), and explains the higher abundance in acidic stands. The Reyranella genus was enriched in lime mono stands, which was consistent with previous reports of forest soil bacterial communities (Felske et al., 1998; Kim et al., 2013).

Betaproteobacteria were represented by Burkholderiales (unidentified Comamonclass="Disease">adaceae genus, Variovorax and Rhizobacter) and Nitrosomonclass="Chemical">pan class="Disease">adales (unidentified Nitrosomonadaceae genus and Nitrosospira) (Supplementary Figures S4, S5; Supplementary Table S3). At total and active community level, Burkholderiales and Nitrosomonadales were more abundant in lime and hornbeam mono stands (6.8 and 9.5%, respectively) compared to beech and oak mono stands (3.8 and 4.5%, respectively). A similar composition was also found in the mixed stands BHL and HOL compared to mixed stands BHO and BOL. Both orders form part of the nitrogen-fixing bacterial community in forests soils and participate in symbiotic relationships with plants (Cherobaeva et al., 2011; Tkacz and Poole, 2015).

Myxococcales and Declass="Chemical">sulfurellales orders dominated within Deltaclass="Chemical">proteobacteria (Figure 3). class="Chemical">pan class="Species">Haliangium was the most abundant genus within Myxococcales at total community level (1.3%) and active community level (6%), followed by Sorangium. Both genera are commonly found in soils, but only members of Sorangium have been isolated from soil (Dawid, 2000; Fudou et al., 2002). Members of both genera exhibit a capacity for producing secondary metabolites with potential pharmaceutical use (Fudou et al., 2001; Li et al., 2014). The Desulfurellales consisted mainly of Desulfurellaceae family members. Desulfurellaceae were more abundant in lime and hornbeam mono stands compared to beech or oak mono stands at total and active community level (Supplementary Table S3). Desulfurellaceae are obligate sulfur-metabolizing thermophiles that contribute to the sulfur cycle (Flores et al., 2012; Wang et al., 2016). Their presence in forest soils is not widely reported but one study has reported Desulfurellaceae in farm soils (Wang et al., 2016).

The Gammaproteobacteria consisted primarily of Xanthomonclass="Disease">adales (class="Chemical">pan class="Species">Rhodanobacter, Acidibacter, and an unidentified genus). Members of Xanthomonadales, including Rhodanobacter, have been reported to prefer soil environments with low pH and high C/N ratios, which promote efficient denitrification (van den Heuvel et al., 2010; Green et al., 2012; Prakash et al., 2012). Correspondingly, Xanthomonadales were enriched at total and active community in beech (6.9 and 5.1%, respectively) and oak mono stands (7.7 and 7.3%, respectively) compared to lime mono stands (5.2 and 3%, respectively).

The composition of Acidobacteria varied considerably across forest stands, through several unidentified subgroups (Supplementary Figures S4, S5). Acidobacteriales, Solibacterales, and subgroup 2 were more abundant in beech and oak mono stands compared to class="Species">lime and hornbeam mono stands (Suclass="Chemical">pclass="Chemical">plementary Figures S4, S5; Suclass="Chemical">pclass="Chemical">plementary Tables S3 and S4). In contrast, subgrouclass="Chemical">p 6 showed higher relative abundance in class="Chemical">pan class="Species">lime and hornbeam mono stands (19% at DNA level and 7.3% at RNA level) than in beech and hornbeam mono stands (6.3% at DNA level and 2.2% at RNA level). In mixed stands, Acidobacteriales, Solibacterales, and subgroup 2 exhibited similar trends as in mono stands and were more abundant in beech and oak mixed stands, BHO and BOL. Subgroup 6 showed higher abundance in lime and hornbeam mixed stands BHL and HOL, at total and active community level. Acidobacteria correlate negatively with pH and nutrient availability (Sait et al., 2006; Clivot et al., 2012; Miyashita, 2015; Kielak et al., 2016a). Acidobacterial groups with copiotrophic lifestyles, which includes subgroup 6, were enriched in more neutral, nutrient-rich environments (Naether et al., 2012; Huber et al., 2016), which confirms our results.

Actinobacteria, which includes plant pathogens and members capable of producing secondary metabolites (Ventura et al., 2007; Barka et al., 2016), was mainly represented by the Gaiellales, Frankiales, Acidimicrobiales, and Solirubrobacterales orders (Supplementary Tables S3 and S4). Frankiales abundance increased in the active community (4%) compared to the total community (2.9%). Many dominant genera within these orders were unclassified but we identified class="Species">Acidothermus within Frankiales. class="Chemical">pan class="Species">Acidothermus showed lower abundance in lime mono stands compared to beech, hornbeam, and oak, at total and potentially active community levels (Supplementary Tables S3 and S4). Acidothermus was previously recorded in high relative abundance (11.5% at DNA level) from a natural forest (Kim et al., 2015). Additionally, the only cultured species, Acidothermus cellulolyticus, was isolated from an acidic hot spring, which points to an acidophilic lifestyle, as shown by our study (Mohagheghi et al., 1986; Barabote et al., 2009).

Representatives of Bacteroidetes belonged predominantly to Cytophagaceae (Cytophagales), Flavobacteriaceae (Flavobacteriales), and Chitinophagaceae (Sphingobacteriales) families, which showed higher abundance in class="Species">lime and hornbeam mono stands than in beech and oak mono stands (Suclass="Chemical">pclass="Chemical">plementary Figures S4, S5; Suclass="Chemical">pclass="Chemical">plementary Tables S3, S4). The genus Flavobacterium within Flavobacteriaceae, is a common soil inhabitant, and was more enriched at active community level. Flavobacterium was more abundant in class="Chemical">pan class="Species">lime and hornbeam mono stands (2.8% at DNA level and 2.9% at RNA level) than in beech and oak mono stands (0.78% at DNA level and 0.56% at RNA level). Species of Flavobacterium have been reported in plant root associations and deadwood, and harbor broad physiological capabilities including lignolytic activity (Dilly et al., 2000; Kolton et al., 2013; Deshmukh and Sao, 2015; Hoppe et al., 2015).

Taxa-Habitat Association Patterns

Bipartite association networks provided insight into bacterial taxa that potentially drive the observed community structures across tree stands (Figure 4; Supplementary Figure S6). Networks were constructed from OTUs showing significant positive associations to specific stands or a combination of stands. We detected 466 and 348 indicator OTUs at total and active community level, respectively. The correlation-based network strongly mirrors the pattern of forest stands recorded during NMDS analysis (Figure 2). Shorter cross edges between beech and oak stands compared to hornbeam or class="Species">lime stands indicate that communities in beech and oak stands are more closely associated with tree sclass="Chemical">pecies than those associated with hornbeam or class="Chemical">pan class="Species">lime stands. No significant positive associations were observed between bacterial genera in oak and lime mono stands.

Figure 4

Association networks between soil bacterial communities (genus level) and mono stands. Source nodes (rounded squares) represent mono species tree stands and edges represent associations between stands and bacterial OTUs (circles, target nodes). Edges are colored according to the source tree species and the length of edges is weighted according to association strength. Unique clusters, which associate with one tree species, consist of nodes colored as the corresponding stand. Numbers of OTUs making up respective unique clusters are given in brackets. Black circles represent OTUs with significant cross association between two or more plots. Target node sizes represent mean relative abundance of OTUs across all mono plots. Data only represents OTUs that showed significant positive association with tree species (p ≤ 0.05). In the case that the genus could not be assigned, the taxonomic name at the highest determined taxonomic resolution is given in parenthesis. For ease of visualization, edges were bundled together, with a stress value of 3. Abbreviations: beech (B), hornbeam (H), lime (L), oak (O), and operational taxonomic unit (OTU).

Association networks between soil bacterial communities (genus level) and mono stands. Source nodes (rounded squares) represent mono species tree stands and edges represent associations between stands and bacterial OTUs (circles, target nodes). Edges are colored according to the source tree species and the length of edges is weighted according to association strength. Unique clusters, which associate with one tree species, consist of nodes colored as the corresponding stand. Numbers of OTUs making up respective unique clusters are given in brackets. Black circles represent OTUs with significant cross association between two or more plots. Target node sizes represent mean relative abundance of OTUs across all mono plots. Data only represents OTUs that showed significant positive association with tree species (p ≤ 0.05). In the case that the genus could not be assigned, the taxonomic name at the highest determined taxonomic resolution is given in parenthesis. For ease of visualization, edges were bundled together, with a stress value of 3. Abbreviations: beech (B), hornbeam (H), pan class="Species">lime (L), oak (O), and oclass="Chemical">perational taxonomic unit (OTU).

Unique clusters representing OTUs associated significantly with only one mono stand, accounted for 58 and 64% of all network OTUs at total and active community level, respectively. The higher number of genera that define unique clusters in class="Species">lime stands in the total (191) and active (169) community arises from the high bacterial diversity observed in class="Chemical">pan class="Species">lime stands compared to hornbeam, beech, and oak stands. Unique clusters in mixed stand were less common and instead, OTUs formed cross associations with two or more mixed stands (Supplementary Figure S6). This provides further evidence that bacterial communities in mixed stands are composed of members associated with corresponding mono stands and illustrates the importance of tree species identity over tree species richness in shaping soil bacterial community.

Genera in unique clusters belonged to dominant phyla and orders described in this work. At total community level, we identified unclassified OTUs from Acidobacteria subgroup 6, Nitrosomonclass="Disease">adacea and class="Chemical">pan class="Disease">Reyranella genus in lime mono stands. Any overlap of significant bacterial cross-associations between beech and oak, and hornbeam and lime mono stands was provided by a few genera and comprise, among others, Rhodanobacter, Rhizomicrobium, Acidothermus, Bryobacter, Granulicella, members of Gemmatimonadacea, and members of Cytophagaceae. This pattern is similarly reflected in the active community but with fewer genera. Analyses of soil environments have revealed shared taxonomic groups, but only a few genera participate in distinguishing one soil habitat from another (Hartmann et al., 2015; Rime et al., 2016). A study on microbial community conversion between organic and conventional farming showed that only 12% of bacterial OTUs constituted the management-specific community, with 49% showing significant unique association to a specific management type (Hartmann et al., 2015). Our results follow this trend, as only a small fraction of total OTUs define the difference between soil bacterial communities of different tree stands.

Bacterial Functional Profiles Across Forest Stands

The 16S rRNA transcript-based communities were used to predict active metabolic processes in forest soil. It has to be noted that rRNA abundance is only a qualitative index for activity but not direct measure of activity (Blazewicz et al., 2013). Thus, the predicted functional profile does not necessarily reflect the direct activity of the studied organisclass="Disease">ms. class="Chemical">pan class="Disease">Additionally, OTUs derived from unknown taxa limit functional predictions. Nevertheless, it has been shown for bacterioplankton and soil bacterial communities, as well as for communities in other environments that 16S rRNA-derived functional prediction are in good agreement with those derived from direct sequencing of corresponding metagenomes and metatranscriptomes (Aßhauer et al., 2015; Kaiser et al., 2016; Wemheuer et al., 2017).

In general, predicted metabolic functional profiles did not follow the stand-specific trend recorded for the bacterial community structures but showed strong grouping according to season (p = 0.001) (Figure 5; Supplementary Table S5). Only functional genes related to class="Chemical">methane metabolism showed association with tree sclass="Chemical">pecies (class="Chemical">p = 0.033) (Suclass="Chemical">pclass="Chemical">plementary Figure S7). This is in accordance with a reclass="Chemical">port showing that methylotroclass="Chemical">phic bacteria correlated with shifts in soil class="Chemical">pH in a beech-dominated deciduous forest (Morawe et al., 2017). In general, gene functions associated with class="Chemical">pan class="Chemical">carbon and nitrogen metabolism in spring were distinct from summer and autumn, which grouped more closely (Figure 5). As a subset of carbon metabolism, gene functions associated with methane metabolism also followed this trend. Gene function for sulfur metabolism also showed separation between seasons. As light intensity and temperature increase in spring (April), trees begin to form leaves and increase photosynthetic productivity, which peaks in summer (July) (Goodale et al., 2015). In autumn (October), light intensity and temperature decline and trees, generally, respond through increased litterfall (Goodale et al., 2015; Žifčáková et al., 2016). However, most litter decomposition takes place during the summer (Sohng et al., 2014) and may continue into autumn. Correspondingly, we concluded that belowground metabolic processes from summer extended into autumn, giving rise to similar bacterial functional profiles.

Figure 5

Non-metric multidimensional scaling (NMDS) of functional genes in key biogeochemical processes. Ordination is based on Bray-Curtis distance matrices of functional genes (KEGG orthologs) predicted by Tax4Fun. A subset of genes predicted for each metabolic process were used to generate distance matrices, which were subsequently summarized by tree species (n = 6).

We observed several predicted functional gene profiles related to class="Chemical">carbon (including class="Chemical">pan class="Chemical">methane) and nitrogen metabolism, key pathways in microbial biogeochemical cycling (Figure 6; Supplementary Table S6). The higher abundance of genes observed for carbon metabolism agree with other studies showing that C cycling is a central part to bacterial metabolism in forest soils (Nacke et al., 2014; Siles and Margesin, 2017). Carbon turnover in soils is mediated by carbohydrate-active enzymes (CAZymes), which act on labile C compounds, recalcitrant cellulose or hemicellulose and fungal biomass (López-Mondéjar et al., 2016; Žifčáková et al., 2017). In our study, gene function related to C fixation, including ribulose 1,5-bisphosphate carboxylase/oxygenase (RuBisCo), showed significantly higher abundance in spring and autumn (p = 0.013). However, the relative abundances of C degradation genes (cellulases, hemicellulases, and chitinases) were higher in summer and autumn than in spring. Due to temperature dependence, CAZymes activity was demonstrated to correlate with seasonal changes, which subsequently impacts C cycling (Žifčáková et al., 2017). The quality of soil C input from root exudates and litter varies seasonally (Siles and Margesin, 2017) and explains the seasonal changes in the relative abundance of C fixation and degradation gene functions. A similar increase in C degradation gene function in autumn compared to spring was also reported in a mixed deciduous and coniferous forest ecosystem (Siles and Margesin, 2017).

Figure 6

Overview of predicted bacterial functional categories involved in key energy pathways. Genes were predicted from KEGG orthologs with Tax4Fun (Aßhauer et al., 2015). Each box plot represents the mean relative abundance of predicted genes for the eight mono and mixed stands (each with six replicates) calculated for spring, summer, and autumn. Genes used are listed in Supplementary Table S4. Mean values with identical letters, determined by ANOVA with Tukey HSD post hoc test, share significant similarity among seasons (p ≤ 0.05).

Overview of predicted bacterial functional categories involved in key energy pathways. Genes were predicted from KEGG orthologs with Tax4Fun (Aßhauer et al., 2015). Each box plot represents the mean relative abundance of predicted genes for the eight mono and mixed stands (each with six replicates) calculated for spring, summer, and autumn. Genes used are listed in Supplementary Table S4. Mean values with identical letters, determined by ANOVA with Tukey pan class="Disease">HSD class="Chemical">post hoc test, share significant similarity among seasons (class="Chemical">p ≤ 0.05).

The ability to utilize a wide range of C sources to obtain metabolic energy is shared among several bacterial taxa across phyla and enables several groups to inhabit the same envclass="Chemical">ironment with different nutrient niches (Llclass="Chemical">pan class="Disease">adó et al., 2017). Bacterial phyla comprise metabolically versatile genera that carry out the same general metabolic processes and are functionally redundant (Burke et al., 2011). Therefore, the relative abundance of genes in an environment does not always correlate with metabolic activity of the corresponding pathways, as some bacteria are able to use more efficient pathways to metabolize the same substrate (Rocca et al., 2015). As C allocation is mediated by bacteria through decomposition of organic matter, particularly through lignocellulose breakdown (Lladó et al., 2017), we also found potential cellulolytic genera in our study. These included Burkholderia, Variovorax, and Flavobacterium (Haichar et al., 2007; Ulrich et al., 2008; Schellenberger et al., 2010). Additional genera with potential cellulolytic potential were classified as rare (less than 1%) in our dataset, and included actinobacterial genera (Arthrobacter, Cellulomonas, Kitasatospora, Oerskovia, Micromonospora, and Streptomyces), alphabroteobacterial genera (Mesorhizobium, Methylobacterium, Sphingomonas), the gammaproteobacterial genus Dyella and genera from Firmicutes (Bacillus, Paenibacillus).

Tax4Fun predicted gene functions encoding class="Chemical">methane monoclass="Chemical">pan class="Chemical">oxygenase (MMO), a key enzyme in the oxidation of methane to methanol, harbored by methylotrophic and methanotrophic bacteria (Supplementary Table S6; Hakemian and Rosenzweig, 2007; Hoppe et al., 2015). Methylotrophs metabolize single carbon substrates, as a by-product of lignin degradation (Hoppe et al., 2015). Potential methylotrophs were observed in Rhizobiales (Methylocella, and Methyloferula), gammaproteobacterial Methylococcales (Methylococcaceae), and Verrucomicrobia (Methylacidiphilum). Seasonal inputs of organic matter from decomposition explain the strong effect of season on the abundance of methanotrophs. For example, Hoppe et al. (2015) demonstrated an increase of methylotrophic Rhizobiales at different stages of deadwood decomposition in a beech and spruce deciduous forest.

The presence of gene functions for assimilatory (sir, class="Chemical">cys genes) and dissimilatory (dsr) class="Chemical">pan class="Chemical">sulphate reduction revealed the potential presence of sulfur metabolizing organisms. However, the relative abundance of these genes did not exhibit significant differences across seasons (p = 0.818; Figure 6). This is explained by reports that sulphate reducing bacteria tend to be rare and are likely CRTs, within the environment despite carrying out the major part of sulfur metabolism (Yousuf et al., 2014b; Hausmann et al., 2016). Candidates detected included Desulfurellales, Desulfurbacterales, and Desulfuromonadales within Deltaproteobacteria, Thiohalophilus, and Thioalkalispira from Gammaproteobacteria, Comamonadaceae (Betaproteobacteria), Rhodopseudomonas (Alphaproteobacteria), and Desulfosporosinus (Firmicutes). Several studies have identified these taxa as major components of sulfur cycling in soil and sediment ecosystems (Baker et al., 2015; Balk et al., 2015; Ling et al., 2015; Hausmann et al., 2016).

Gene functions such as nir, nif, hao, and amo for the main processes in class="Chemical">nitrogen cycling (nitrification, class="Chemical">pan class="Chemical">nitrate reduction, and anammox) were observed (Supplementary Table S6; Norton et al., 2002; Giles et al., 2012). Groups associated with N cycling from our study potentially belong to Nitrospira, all genera within Nitrosomonadaceae and Rhizobiales (Bradyrhizobium and Rhizobium), Rhodospirillaceae (Azospirillum), and Actinobacteria (Arthrobacter). Gene function for N metabolism was significantly more abundant during spring and autumn (p ≤ 0.01) (Figure 6), and similar to C cycling genes, corresponding to periods of increased nutrient availability. Taxonomic analysis revealed that the diversity of N cycling bacterial community members spans several taxonomic groups. This also demonstrates the complexity of bacteria-mediated N cycling, which is closely interlinked with C cycling (Yousuf et al., 2014a). Thus, bacteria have developed alternative pathways enabling facultative metabolism of both carbon and nitrogen substrates (Freedman et al., 2013; Yousuf et al., 2014a; Rocca et al., 2015). Furthermore, the shared trend in seasonal abundance observed for C and N metabolic gene functions in our study was expected, as soil C allocation affects soil pH and C/N ratio, which subsequently impact nitrification, denitrification, and C cycling (Townsend et al., 2011; Cardenas et al., 2018).

Conclusions

Tree species exhibited a strong positive correlations with soil bacterial diversity and composition at entire and active soil bacterial community level, which supported our first hypothesis (1) that tree species identity drives bacterial community structure at entire and potentially active bacterial community level. This potentially is a result of both direct and indirect factors such as litterfall and root exudates, which change soil pH, C/N ratio, N and P availability. Beech and oak mono stands displayed low pH and high C/N ratio and correspondingly, showed a higher relative abundance of oligotrophic and lower relative abundance of copiotrophic bacterial taxa compared to class="Species">lime and hornbeam mono stands. Mono stands comclass="Chemical">pared to mixed stands showed a higher number of indicator OTUs corresclass="Chemical">ponding to organisclass="Chemical">pan class="Disease">ms closely associated with each stand. Most indicator OTUs belonged to Rhizobiales, indicating the widespread physiological adaptation of its members to different environments. Our second hypothesis that (2) metabolic functions are also driven by tree species was not supported. We did not observe stand-specific effects on predicted bacterial metabolic functions, except for gene functions related to methane metabolism. Interestingly, predicted functional metabolic profiles correlated significantly with season. This was attributed to functional redundancy across different taxonomic groups. We observed more gene functions associated with C fixation and degradation (including methane metabolism), compared to nitrogen metabolism. This supports evidence that forest ecosystems play a central role in carbon storage and contribute to global carbon cycling. The abundance of C and cycling genes showed similar increase in spring and autumn, which was linked to shared metabolic pathways across different bacterial taxa. An increase in gene functions for both processes corresponded to periods of increased soil nutrient availability, as a response to increased root productivity (spring) and litterfall (autumn). We demonstrated that hypothesis (3) soil physicochemical parameters are strongly influenced by tree species identity (mono stands) and richness (mixed stands), which subsequently drive the observed differences in bacterial community structure at total and active community level. Furthermore, our hypothesis (2) that soil bacterial structure across forest stands is shaped by season to a lesser extent than by tree species was not fully supported. Season showed a non-significant effect on soil bacterial community composition at total and active community level. The approach to evaluate soil bacteria at total and active community level provided a comprehensive overview of compositional and potential functional changes in soil bacterial communities. Moreover, the ability to monitor taxonomic and functional relationships in individual microbial taxa provides insights into specific impacts of trees on shaping soil microbial communities and improves our understanding of how potential conversion of forest stands effects change in soil microbial community and function.

Author Contributions

RD designed and conceived the study. Soil sampling for prokaryotic community analysis was performed by class="Disease">AD. class="Chemical">pan class="Disease">AD carried out the field and laboratory work. MS and EV set up experimental design and site and provided data. AD and DS prepared and analyzed the data. AD generated the draft version of the manuscript. All authors interpreted the results and contributed to writing the final version of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.