Soil multifunctionality is affected by the soil environment and by microbial community composition and diversity.

Microorganisms are critical in mediating carbon (C) and nitrogen (N) cycling processes in soils. Yet, it has long been debated whether the processes underlying biogeochemical cycles are affected by the composition and diversity of the soil microbial community or not. The composition and diversity of soil microbial communities can be influenced by various environmental factors, which in turn are known to impact biogeochemical processes. The objectives of this study were to test effects of multiple edaphic drivers individually and represented as the multivariate soil environment interacting with microbial community composition and diversity, and concomitantly on multiple soil functions (i.e. soil enzyme activities, soil C and N processes). We employed high-throughput sequencing (Illumina MiSeq) to analyze bacterial/archaeal and fungal community composition by targeting the 16S rRNA gene and the ITS1 region of soils collected from three land uses (cropland, grassland and forest) deriving from two bedrock forms (silicate and limestone). Based on this data set we explored single and combined effects of edaphic variables on soil microbial community structure and diversity, as well as on soil enzyme activities and several soil C and N processes. We found that both bacterial/archaeal and fungal communities were shaped by the same edaphic factors, with most single edaphic variables and the combined soil environment representation exerting stronger effects on bacterial/archaeal communities than on fungal communities, as demonstrated by (partial) Mantel tests. We also found similar edaphic controls on the bacterial/archaeal/fungal richness and diversity. Soil C processes were only directly affected by the soil environment but not affected by microbial community composition. In contrast, soil N processes were significantly related to bacterial/archaeal community composition and bacterial/archaeal/fungal richness/diversity but not directly affected by the soil environment. This indicates direct control of the soil environment on soil C processes and indirect control of the soil environment on soil N processes by structuring the microbial communities. The study further highlights the importance of edaphic drivers and microbial communities (i.e. composition and diversity) on important soil C and N processes.

Introduction

Soils harbor an enormous diversity of microorganisms, among which bacteria, archaea and fungi play pivotal roles for ecosystem functioning, such as regulating organic matter decomposition and soil C dynamics, and mediating nutrient cycling (Bardgett et al., 2008; Singh et al., 2010; Wagg et al., 2014). As microbial habitats, the soil environment was reported to exert substantial impacts on microbial community structure and diversity (Lauber et al., 2008; Rasche et al., 2011; Richter et al., 2018). A broad range of edaphic variables such as soil pH, texture, moisture, temperature, organic C and nutrient content were recognized to influence the composition and diversity of soil microbial communities (Brockett et al., 2012; Cookson et al., 2007; Rousk et al., 2010a). At the global scale, soil pH is regarded as the key predictor of soil bacterial community composition and diversity (Fierer and Jackson, 2006; Rousk et al., 2010a; Zhalnina et al., 2015). Soil texture, particularly clay and silt content, is closely related to soil organic C (SOC) content and nutrient availability, and was shown as another key driver of microbial community composition and diversity (Hansel et al., 2008; Kallenbach et al., 2016). Microbial community composition and diversity could also be affected by soil nutrient content (Koyama et al., 2014; Pan et al., 2014); N and P addition for example were reported to increase bacterial to fungal phospholipid fatty acid ratios (Dong et al., 2015) and change microbial diversity (Leff et al., 2015; Ling et al., 2017). Beyond the recorded influence of the soil environment on microbial communities, there is a wealth of studies on edaphic and environmental effects on soil functions such as soil formation, organic matter decomposition and substrate use efficiency (e.g. Bonner et al., 2018; Borken and Matzner, 2009; Colman and Schimel, 2013; Davidson et al., 1998; Davidson and Janssens, 2006; Hu et al., 2018). For instance, temperature, soil moisture, substrate availability and nutrient limitations were suggested to affect soil C metabolism including microbial growth, respiration, C use efficiency and microbial biomass turnover (Dijkstra et al., 2015; Hagerty et al., 2014; Manzoni et al., 2012; Schindlbacher et al., 2015; Takriti et al., 2018; Zheng et al., 2019). Soil organic nitrogen (N) transformations such as gross protein depolymerization, gross N mineralization and gross nitrification rates can be controlled by temperature, soil pH, resource or enzyme availability and substrate quality (Booth et al., 2005; Cookson et al., 2007; Noll et al., 2019; Rustad et al., 2001; Wallenstein and Weintraub, 2008; Wanek et al., 2010). Soil functions can also be driven by soil microbial community composition and diversity (Balser and Firestone, 2005; Bonner et al., 2018; Creamer et al., 2015; Don et al., 2017; Schimel and Schaeffer, 2012). Microbial growth and CUE were found to be influenced by bacterial versus fungal dominance (Soares and Rousk, 2019). Soil ammonia-oxidizer populations such as bacterial and archaeal nitrifiers can promote gross nitrification rate (Li et al., 2018; Prommer et al., 2014; Stieglmeier et al., 2014). Soil enzymes activities were found to be shaped by microbial communities (Gallo et al., 2004; Schnecker et al., 2015; Waldrop et al., 2000). It is undoubtedly important to study links between microbial communities and single soil functions, which provide valuable information on microbial drivers of specific processes.

Despite considerable research efforts made into examining how microbial community composition and diversity drives single soil functions, in recent years there is an emerging field of research began to investigate how microbial communities maintain ecosystem multifunctionality based on both observational and manipulative studies (Bastida et al., 2016; Delgado-Baquerizo et al., 2017a, 2017b; 2016; Wagg et al., 2014). These studies calculated multifunctionality indices and attempted to investigate how microbial community composition, richness or diversity drive such multifunctionality. Some also accounted for environmental factors (Delgado-Baquerizo et al., 2016; Thakur et al., 2018), but most of them focused on dissecting the individual effects of single edaphic factors, and few of them have regarded edaphic factors as an integral construct to represent the multivariate soil environment, nor investigated combined effects of the soil environment on microbial community composition, diversity and C and N processes together. Moreover, most studies focused only on effects of bacterial communities and fewer considered effects of fungal or archaeal communities (alongside bacterial ones) on various soil processes, and thus the latter effects remain elusive, specifically across different soils and for multiple processes (Graham et al., 2016). Hence few studies have investigated the soil environment, microbial community structure, richness and diversity, extracellular enzyme patterns and soil C and N processes cohesively. This represents a major knowledge gap given that soils are complex systems that encompass a wide variety of abiotic and biotic characteristics, which means that no single soil parameter can explain single or multiple soil processes alone. It is therefore important not only to assess the influence of single edaphic factors but also of the combined effects of multiple factors on multiple processes to allow firm conclusions on the environment-microbial community-function coupling. This is particularly important as soils only provide their ecosystem services based on their multifunctional integrity (Delgado-Baquerizo et al., 2016; Wagg et al., 2014).

The objectives of this study were to test the effects of multiple edaphic drivers, microbial community composition and diversity on soil multiple functions (i.e. soil enzyme activities, soil C and N processes). Towards this end, we examined bacterial/archaeal and fungal community composition using DNA-based sequencing methods (Illumina MiSeq) and linked them to a series of edaphic variables, as well as a wide range of soil processes and soil enzyme activities. We studied individual effects of single edaphic factors on microbial community composition, soil process rates and extracellular enzyme activities, and also investigated the combined effects of soil parameters on soil multiple processes as matrices by Mantel tests, which provides a more comprehensive understanding of using environmental and microbial data to predict the multiple functions of soil ecosystems. Soils from three land uses (cropland, grassland and forest) deriving from two bedrock forms (silicate and limestone) were collected to test for the generality of the patterns.

Materials and methods

Site description and soil sampling

Soils were collected from three land-use areas (cropland, grassland and forest) at two sites in the central Enns valley, Styria, Austria: (1) LFZ Raumberg-Gumpenstein (47° 29′ N, 14° 6′ E, 690m a.s.l) and (2) Moarhof in Trautenfels-Pürgg (47° 30′ N, 14° 4′ E, 708m a.s.l) in June 2016. These sites were located on opposite sides of a valley, thus experiencing similar climate, with mean annual precipitation of 980mm and mean annual temperature of 7.2 °C. They differed in their associated bedrock material: soils from LFZ Raumberg-Gumpenstein derived from silicate bedrock (Gneiss), whereas Moarhof soils derived from calcareous bedrock (limestone, dolomite). Overall soils were classified as Luvisols (limestone, L) and Cambisols (silicate, S) (SWSR, 2015).

Croplands (C) were cultivated with a mixture of vegetables including cabbage, bean, potato, and onion on silicate soils and by a mixture of oat, barley and wheat on limestone soils. Grasslands (G) were permanent grasslands grazed by sheep on silicate sites and by cattle on limestone sites. Forests (F) were dominated by spruce (Picea abies L.) on silicate sites and by spruce and ash (Fraxinus excelsior L.) on limestone sites. After removing the litter and organic layers, four independent replicates of mineral soils were sampled to a soil depth of 15 cm using a root corer (Eijkelkamp, Netherlands) with 7.5 cm in diameter at each site. The four biological replicates were sampled to cover major heterogeneities in site topography (bottom of slope, upper and lower slope, hilltop) or crop plants. All soils were sieved to 2mm and stored at 4 °C until further analysis. Soil replicates were stored and processed independently.

Soil physicochemical and biological analysis

Soil texture, bulk density, cation exchange capacity (CEC), base saturation, carbonate content, exchangeable Ca2+, Mg2+, K+, Na+, Al3+, Fe3+, Mn2+ and H+ were determined by the soil analysis laboratory of the Federal Office for Food Safety (AGES, Vienna, Austria) according to standard protocols. Aliquots (10 g) of fresh soils were dried in a drying oven at 80 °C for three days to determine the soil water content (SWC). Soil pH was measured in Milli-Q water (soil: solution ratio = 1:2.5 (w: v)) using an ISFET electrode (Sentron, Netherlands). Total soil organic C (SOC) and total N (TN) were analyzed in aliquots of oven-dried and ball milled (MM200, Retsch, Germany) soils using an Elemental Analyzer (Carlo Erba 1110, CE Instruments) coupled to a DeltaPlus Isotope Ratio Mass Spectrometer (Finnigan MAT, Thermo Fisher, Germany) via a Conflo III interface (Thermo Fisher, Austria), after removing carbonate using 2M HCl. Soil dissolved organic C (DOC) and total dissolved N (TDN) were measured by a TOC/TN analyzer (TOC-VCPH/TNM-1, Shimadzu, Austria) in 1M KCl (1:7.5 (w: v) for 60 min) extracts. Ammonium (NH4+) and nitrate (NO3−) concentrations were determined in the same extracts photometrically (Hood-Nowotny et al., 2010). Soil total P (TP) and total inorganic P (TIP) were measured in 0.5M H2SO4 extracts of ignited (450 °C, 4 h) and control soils (Kuo, 1996) by malachite green measurements of reactive phosphate (Lajtha et al., 1999). Soil total organic P (TOP) was calculated as the difference between TP and TIP. Dissolved inorganic P (DIP) was determined using malachite green in 0.5M NaHCO3 (pH 8.5; 1:7.5 (w: v)) extracts after acidification with H2SO4. Acid persulfate digestion (Lajtha et al., 1999) was applied to measure total dissolved P (TDP) and allowed calculating dissolved organic P (DOP). Soil microbial biomass C (MBC) and microbial biomass N (MBN) were determined using chloroform fumigation extraction (Vance et al., 1987) for 48 h. Soil microbial biomass P (MBP) was also performed by chloroform-fumigation extraction but using 0.5M NaHCO3 instead of 1M KCl extractions (Brookes et al., 1982), and was calculated as the difference in TDP in extracts of fumigated and non-fumigated soils.

Determination of soil process rates and extracellular enzyme activities

All soil process rates and extracellular enzyme activities were measured after the soils were pre-incubated at 20 °C and 60% water holding capacity (WHC) for 7 days, and the values were normalized to MBC as the sequencing data of bacterial/archaeal and fungal communities only provide relative and not absolute abundances of micro-organisms. Soil microbial growth, basal respiration and C use efficiency (CUE) were determined based on an 18O-water method (Zheng et al., 2019). N use efficiency (NUE) (Mooshammer et al., 2014) and gross rates of N transformation processes including gross protein depolymerization, gross N mineralization and gross nitrification rates were determined using isotope pool dilution (IPD) assays (Wanek et al., 2010).

Microplate fluorimetric and photometric assays were applied to measure potential extracellular enzyme activities in soils (Kaiser et al., 2010). Phenoloxidase (POX) activities were measured photometrically using ABTS (2,2′-azinobis-(-3 ethylbenzothiazoline-6-sulfonic acid)) as the substrate (Zheng et al., 2019; Floch et al., 2007). β-glucosidase (BG) and phosphatase activities were measured photometrically using p-nitrophenyl (pNP)-linked β-glucopyranoside and p-nitrophenyl (pNP)-linked phosphate as the substrates (Robertson et al., 1999), respectively. L-leucine-7-amido-4-methyl coumarin (AMC) was used as a substrate to measure leucine amino peptidase activities fluorimetrically (Hu et al., 2018).

DNA extraction, amplification and sequencing

The four soil biological replicates that were pre-incubated individually at 20 °C and 60% WHC for 7 days were frozen at −80 °C for further DNA analysis. Soil DNA was extracted from 0.4 g of frozen soil using the FastDNA™ SPIN Kit (MP Biomedicals, Germany) following the modified manufacturers’ recommendations (Spohn et al., 2016). Aliquots (50 μl) of DNA extracts were purified with OneStep™ PCR Inhibitor Removal Kit (Epigenetics, USA) and subsequently quantified using a microtiter plate assay with Quant-iT™ PicoGreen® dsDNA Reagent (Thermo Fisher, Germany). DNA templates of each sample were prepared by diluting the purified DNA to 10 ng μl−1 with nuclease-free water (Carl Roth, Germany).

The V4 hypervariable region of the 16S rRNA gene and the internal transcribed spacer 1 (ITS1) region were amplified via polymerase chain reaction (PCR) using multiplexed barcoded amplicon sequencing on the Illumina MiSeq platform (Illumina, San Diego, CA, USA), as described by Herbold et al. (2015). Briefly, the V4 hypervariable region of the 16S rRNA gene (ca. 292 bp) was amplified using the modified 515-F (5′ GTG CCA GCM GCC GCG GTA A 3′) and modified 806-R (5′ GGA CTA CHV GGG TWT CTA AT 3′) primer pairs (Caporaso et al., 2012), while the ITS1 region (ca. 350 bp) was amplified using the ITS1F (5′ CTT GGT CAT TTA GAG GAA GTA A 3′) and ITS2 (5′ GCT GCG TTC TTC ATC GAT GC 3’) primer pairs (Gardes and Bruns, 1993; White et al., 1990). A second PCR reaction with primers containing sample-specific barcodes was used that allowed pooling the libraries of target genes (Herbold et al., 2015). PCR amplification of each soil replicate was performed in triplicate in 20 μL of cocktail, which consisted of 14.32 μL nuclease-free water, 2 μL 1 x Dream Taq Green Buffer, 2 μL 0.2mM of nucleotide dNTP mixture, 0.08 μL 0.08 μg μl−1 of BSA, 0.5 μL 0.25 μM of each primer, 0.1 μL 1.25 U of DreamTaq Green DNA Polymerase (all from Thermo Fisher Scientific, Waltham, United States) and 1 μL of DNA template (10 ng per reaction). Following amplification, the triplicate of first-step PCR products for each soil replicate were pooled and purified using the ZR-96 DNA Clean-up kit™ (Zymo, United States) and 3 μL of the purified samples were used for the second PCR step (50 μL in total) which contained 32.6 μL nuclease-free water, 5 μL 1 x Dream Taq Green Buffer, 5 μL 0.2mM of nucleotide dNTP mixture, 0.2 μL 0.08 μg μl−1 of BSA, 4 μL 0.8 μM of barcode primer, and 0.25 μL 1.25 U of DreamTaq Green DNA Polymerase. The barcodes can be found in Table S1.

For 16S rRNA gene amplification, the following thermocycling conditions were carried out for the first PCR step: 94 °C for 4 min, followed by 25 cycles of 94 °C for 30 s, 52 °C for 45 s, 72 °C for 45 s, and a final step at 72 °C for 10 min. To amplify the ITS1 region, the thermocycling conditions for the first PCR step were: 94 °C for 3 min, followed by 23 cycles of 94 °C for 45 s, 52 °C for 60 s, 72 °C for 90 s, and a final step at 72 °C for 10 min. PCR products were purified using the ZR-96 DNA Clean-up kit™ (Zymo, United States) and the correct size of PCR products was verified using 1% agarose gel electrophoresis, visualized with GelRed Nucleic Acid Stain (Biotium, United States). The concentration of DNA of samples was quantified with the Quant-iT™ PicoGreen® dsDNA Reagent (Thermo Fisher, Germany) assay and pooled at equal molarity of 20 × 109 copies per sample. Negative controls for both DNA extraction and PCR were performed throughout the experiment. Sequencing was performed by Microsynth AG (Balgach, Switzerland) on a MiSeq platform (Illumina, United States). The library was prepared by adaptor ligation and PCR using the TruSeq Nano DNA Library Prep Kit (Illumina, United States) according to the TruSeq nano protocol (Illumina, FC-121-4003), but excluding the fragmentation step. The MiSeq was run in the 2 × 300 cycle configuration using the MiSeq Reagent kit v3 (Illumina). The sequence data were deposited into the NCBI Short Read Archive under BioProject accession number PRJNA551019.

Sequencing data processing

Paired-end MiSeq reads were demultiplexed for both 16S rRNA and fungal ITSx as described previously (Herbold et al., 2015). Paired end reads were end-trimmed (Q = 10) and further end-trimmed to 200 nt to facilitate merging with fastq-join (Aronesty, 2013). The bacterial/archaeal OTUs (16S rRNA gene) were identified using a 97% identity threshold and OTU representatives were classified using Mothur's implementation of a Naïve Bayesian sequence classifier (Schloss et al., 2009; Wang et al., 2007) and the taxonomic assignment was performed using the SILVA 119 SSU NR99 database (Quast et al., 2013), at a confidence cutoff of 80%.

The paired-end raw fungal ITS1 sequence MiSeq reads were extracted using ITSx prior to clustering data into OTUs (Bengtsson-Palme et al., 2013), singletons were removed and fungal OTUs were identified using a 99% identity threshold. OTUs were clustered, checked for chimeras and assembled into incidence tables using Uparse (Edgar, 2013). The taxonomic assignment for the ITS1 region was implemented by a Naïve Bayesian sequence classifier (Wang et al., 2007) along with the Warcup training set Version 2 (Deshpande et al., 2016) at a confidence cutoff of 80%.

All statistical analyses were performed using R software version 3.4.3 (R Core Team, 2017). The data of edaphic properties, soil process rates and extracellular enzyme activities as well as sequence data were log-transformed when necessary to improve normality and homogeneity. Two-way ANOVAs and Tukey-HSD were applied to test for the effects of bedrock and land use on selected edaphic properties, extracellular enzyme activities and soil biogeochemical process rates. One-way ANOVA was applied to test land use effects on the relative abundance of selected fungal classes (i.e. Agaricomycetes and Eurotiomycetes). Regression analysis among or across selected edaphic parameters, process rates and soil enzyme activities were applied and expressed as Pearson coefficient (R). Spearman correlations were computed between the relative abundances of microbial phyla/classes and edaphic variables and soil process rates. As only a very low diversity and relative abundance of archaea were detected across the tested soils, bacteria/archaea were analyzed together based on sequencing of the 16S rRNA gene. Amplicon libraries were rarefied to 1410 and 744 sequences for bacteria/archaea and fungi, respectively, to ensure even sampling depth for microbial community relative abundance comparison. The beta-diversity analysis was performed using the ‘vegan’ package. Our focus is mainly on microbial community composition (beta-diversity), but we also analyzed microbial richness (i.e. Chao1) and diversity (i.e. Shannon) indices using the ‘phyloseq’ package. The Shannon index takes both microbial richness and evenness into account. In order to obtain the overall variance in microbial composition, the similarities in OTU composition of bacterial/archaeal and of fungal communities across samples were visualized by nonmetric multidimensional scaling (NMDS) ordinations based on Bray-Curtis dissimilarity. Significant environmental variables with P < 0.05 based on permutated data were selected based on all measured soil parameters and fitted onto the NMDS ordination space using the ‘envfit’ function in the ‘vegan’ R package, and significances of correlations were tested with 999 permutations. The first NMDS dimension was rotated parallel to an external edaphic variable (CEC for 16S rRNA gene NMDS plot, TP for ITS NMDS plot). Function ‘bioenv’ in the ‘vegan’ R package was also used to identify the subsets of edaphic variables that best predicted microbial community composition (Clarke and Ainsworth, 1993), and the potential relationship between microbial community composition and the identified subsets of edaphic variables were inferred through canonical correspondence analysis (CCA, Fig. S1 A & B). By using the normalized edaphic data, principal component analysis (PCA) was used to cluster soil samples based on significant edaphic properties that selected by ‘envfit’ function. Analysis of similarities (ANOSIM) with 999 permutations based on Bray-Curtis dissimilarity was conducted to identify differences in OTU composition within and between sites. Similarity percentages of microbial community composition between land uses or bedrocks were calculated using Bray-Curtis dissimilarities with the ‘simper’ function in the ‘vegan’ R package. Mantel tests were performed to evaluate Spearman rank correlations between each two distance matrices or between one single factor and a matrix, and partial Mantel tests were performed to test the Spearman rank correlations between two matrices while controlling for the effect of other matrices. The distances among edaphic variables, soil processes and enzyme patterns were calculated based on Euclidean dissimilarities, while Bray-Curtis distance was used to evaluate dissimilarities among microbial community composition. The selected edaphic variables, extracellular enzyme activities and soil processes were normalized to zero mean and one unit s.d. Finally, all reported P values, except for those calculated by ANOVA and Tukey-HSD tests, were corrected for multiple testing using the ‘p.adjust’ function (method = ‘BH’, n = 798) (Benjamini and Hochberg, 1995) in R.

Results

Soil physicochemical properties, extracellular enzyme activities and soil process rates

We performed two-way ANOVAs and Tukey-HSD to test for the effects of bedrock and land use on soil physicochemical properties, extracellular enzyme activities and soil biogeochemical process rates (Table 1). Silicate soils exhibited lower soil pH, base saturation and CEC as compared to calcareous soils. Bedrock had a weak but significant influence on clay content, while land use strongly impacted it, with highest values in forest soils. TN was much lower in silicate soils compared to calcareous soils. Forest soils had lower TP content than croplands and grasslands. Only phosphatase activities normalized to MBC were significantly higher in silicate soils compared to calcareous soils, while the potential activities of β-glucosidase, phenoloxidase and aminopeptidase, all normalized to MBC, were not influenced by land use or bedrock. Microbial CUE and qCO2 were not affected by land use or bedrock. qGrowth was highest in silicate soils and in forest soils. Microbial NUE, the MBC normalized gross protein depolymerization and gross mineralization rates were all significantly influenced by bedrock, with higher NUE in calcareous soils and higher qGross protein depolymerization and qGross mineralization rates in silicate soils. Gross nitrification rates normalized to MBC were much higher in cropland soils than in grassland or forest soils. Further detailed information on other soil physicochemical and biological properties and soil process rates can be found in Table S2.

Table 1

Edaphic variables, extracellular enzyme activities and soil processes (means ± 1SE, n = 4).

Soil	SC	SG	SF	LC	LG	LF	Two-way ANOVAs			Tukey-HSD
Bedrock	Silicate	Silicate	Silicate	Limestone	Limestone	Limestone	Land use	Bedrock	Land use x bedrock	Cropland	Grassland	Forest	Silicate	Limestone
Land use	Cropland	Grassland	Forest	Cropland	Grassland	Forest
pH (water)	5.90 ± 0.37	5.38 ± 0.18	4.05 ± 0.10	8.15 ± 0.22	6.43 ± 0.19	6.13 ± 0.08	***	***	*	a	b	c	b	a
Base saturation (%)	93.4 ± 3.5	82.7 ± 3.0	5.2 ± 0.2	100.0 ± 0.0	99.6 ± 0.2	99.3 ± 0.1	***	***	***	a	b	c	b	a
CEC (cmol kg−1)	8.38 ± 2.02	5.36 ± 0.25	9.75 ± 0.75	33.7 ± 0.42	22.73 ± 3.23	22.76 ± 0.51	**	***	*	a	b	ab	b	a
Clay (%)	8.50 ± 0.29	5.55 ± 0.24	17.43 ± 0.15	14.80 ± 0.98	8.10 ± 2.17	14.65 ± 1.53	***	*	**	b	c	a	b	a
SOC (mg g−1)	21.8 ± 1.1	26.7 ± 0.9	49.9 ± 7.6	47.0 ± 0.9	47.9 ± 7.6	36.8 ± 2.4	ns	**	***	–	–	–	b	a
TN (mg g−1)	2.18 ± 0.10	2.82 ± 0.11	2.54 ± 0.37	4.77 ± 0.07	4.59 ± 0.58	3.25 ± 0.10	*	***	*	ab	a	b	b	a
TP (mg g−1)	1625 ± 112	1608 ± 56	525 ± 23	1477 ± 31	1555 ± 210	585 ± 84	***	ns	ns	a	a	b	–	–
DOC (μg g−1)	64.6 ± 3.6	85.6 ± 4.3	160.9 ± 13.1	52.5 ± 1.8	53.8 ± 12.6	53.4 ± 6.9	***	***	***	b	b	a	a	b
q β-Glucosidase (nmol (μg MBC)−1h−1)	0.90 ± 0.23	1.06 ± 0.11	1.23 ± 0.51	0.66 ± 0.06	0.59 ± 0.05	0.84 ± 0.03	ns	ns	ns	–	–	–	–	–
qPhenoloxidase (nmol (μg MBC)−1h−1)	9.07 ± 2.64	7.93 ± 0.64	22.37 ± 9.63	1.48 ± 0.16	5.77 ± 1.32	8.74 ± 1.04	ns	ns	ns	–	–	–	a	b
qAminopeptidase (nmol (μg MBC)−1h−1)	0.05 ± 0.02	0.08 ± 0.01	0.05 ± 0	0.10 ± 0.02	0.04 ± 0	0.06 ± 0	ns	ns	**	–	–	–	–	–
qPhosphatase (nmol (μg MBC)−1h−1)	1.30 ± 0.41	2.30 ± 0.28	3.03 ± 0.17	1.05 ± 0.10	1.52 ± 0.24	2.78 ± 0.20	***	ns	ns	c	b	a	–	–
CUE	0.69 ± 0.04	0.71 ± 0.05	0.81 ± 0.06	0.65 ± 0.12	0.74 ± 0.04	0.55 ± 0.07	ns	ns	ns	–	–	–	–	–
qCO2 (ng CO2–C (μg MBC)−1 h−1)	3.79 ± 0.32	3.62 ± 0.35	6.65 ± 0.88	2.26 ± 0.24	2.21 ± 0.09	1.80 ± 0.09	ns	ns	ns	–	–	–	–	–
qGrowth (ng C (μg MBC)−1 h−1)	1.87 ± 0.60	1.45 ± 0.31	1.86 ± 0.98	1.40 ± 0.52	0.77 ± 0.15	1.63 ± 0.50	*	***	***	b	b	a	a	b
NUE	0.72 ± 0.02	0.83 ± 0.03	0.91 ± 0.02	0.98 ± 0	0.96 ± 0	0.79 ± 0.03	ns	***	***	–	–	–	b	a
qGross protein depolymerization (ng N (μg MBC)−1d−1)	534 ± 251	493 ± 73	435 ± 68	51 ± 7.77	86 ± 13	191 ± 12	ns	**	ns	–	–	–	a	b
qGross mineralization (ng N (μg MBC)−1d−1)	5.42 ± 1.77	4.53 ± 0.55	2.96 ± 0.41	0.30 ± 0.02	0.55 ± 0.07	3.20 ± 0.63	ns	***	*	–	–	–	a	b
qGross nitrification (ng N (μg MBC)−1d−1)	2.94 ± 0.67	2.57 ± 0.55	0.68 ± 0.56	3.21 ± 0.90	1.84 ± 0.41	1.86 ± 0.19	*	ns	*	a	a	b	–	–

Significance levels of two way ANOVA: ***, P < 0.001; **, P < 0.01; *, P < 0.05; ns, not significant.

Bacterial/archaeal and fungal community composition and diversity

Sequencing generated a total of 229981 (16S rRNA gene) sequences of 24 sites that were classified into 6989 distinct bacterial/archaeal OTUs, of which 16895 sequences (7.3% of 16S rRNA gene) and 83 OTUs belonged to archaea. After rarefying to identical sequence depth (1384), 31832 high-quality sequences forming 3657 OTUs of bacteria/archaea were retrieved from 23 soil samples (one silicate forest sample (SF1) that yielded less than 1384 sequences was excluded), of which the identified OTUs could be assigned to 32 phyla. A total of 53207 (ITS1 region) sequences generated of 24 sites were classified into 2875 fungal OTUs. After rarefying to identical sequence depth (586), we obtained 13478 high-quality sequences within 1490 OTUs from 23 samples (one silicate forest sample (SF4) that yielded less than 744 sequences was excluded) of the fungal community, of which the identified OTUs belonged to 4 phyla and 19 classes. Microbial richness was estimated by the Chao1 index. The average bacterial/archaeal richness was 1122 ± 83 OTUs and fungal richness was 251 ± 33 OTUs across all soils (Fig. S2 A & B). Microbial diversity was estimated based on the Shannon index. The Shannon index of bacteria/archaea and fungi were 5.79 and 3.84, respectively.

Across the studied soils, the bacterial phylum Proteobacteria (28.9%, consisting of Alpha- (12.7%), Beta- (6.9%), Delta- (5.4%) and Gamma-proteobacteria (3.4%)), Acidobacteria (17.3%), Actinobacteria (10.9%), Bacteroidetes (6.7%) and Verrucomicrobia (5.8%), as well as the archaeal phylum Thaumarchaeota (7.1%) encompassed the largest proportion of sequences of bacterial/archaeal communities (Fig. 1A). Members of the Ascomycota (49.1%), Basidiomycota (22.1%), Zygomycota (1.5%) and Chytridiomycota (1.0%) were prevalent fungal groups across the investigated soils. Of the Basidiomycota phylum, the class Agaricomycetes was the most abundant (21.0%); while of the Ascomycota phylum, the class Eurotiomycetes were the most abundant (17.4%), followed by Sordariomycetes (5.8%), Pezizomycetes (4.6%), Leotiomycetes (3.1%) and Dothideomycetes (2.6%) (Fig. 1B).

Fig. 1

Relative abundance of the bacterial/archaeal (A) and fungal (B) communities at the phylum and the class level respectively based on Illumina sequencing of the 16S rRNA gene and the ITS1 region from three land uses (C: cropland; G: grassland; F: forest) on two bedrock forms (S: silicate; L: limestone). ‘Others’ include phyla or classes with < 1% average relative abundance. Data represent the mean of 4 replicate samples, with the exception of SF (3 replicates).

Edaphic factors shape microbial community composition in soils

Bacterial/archaeal (ANOSIM: R = 0.920, P < 0.001) and fungal (ANOSIM: R = 0.948, P < 0.001) communities clustered significantly by site. Bacterial/archaeal communities were more different between bedrocks (ANOSIM: R = 0.611, P < 0.001) than between land uses (ANOSIM: R = 0.231, P < 0.01). The links between edaphic variables and microbial community composition were visualized in the NMDS ordination plot. The differences in bacterial/archaeal communities were mainly associated with CEC, TN, SOC and pH, with Acidobacteria, Actinobacteria, Thaumarchaeota and Verrucomicrobia accounting for 56% of the overall dissimilarities between silicate and limestone soils (Table S3A). Fungal communities showed stronger differences between land uses (ANOSIM: R = 0.446, P < 0.001) than between bedrocks (ANOSIM: R = 0.255, P < 0.01); forest soils showed higher differences in their fungal communities compared to croplands and grasslands, with Agaricomycetes and Eurotiomycetes explaining 48% of the total dissimilarity between croplands and forests, and 52% between grasslands and forests (Table S3B). The differences in fungal community composition were linked to differences in TP and clay content, with higher TP in cropland and grassland soils compared to forest soils while the latter were high in clay content. Similar edaphic factors were found to be significantly correlated with both bacterial/archaeal and fungal community composition when overlain onto the NMDS ordination space (Fig. 2) including pH, base saturation, CEC, clay, SOC, TN, TP and DOC content, which was corroborated by Mantel tests (Table 2). Based on Mantel tests, bacterial/archaeal communities exhibited significant and strong correlations with soil pH, base saturation, DOC, CEC and TN, followed by relatively weak but significant correlations with SOC, clay and TP. Fungal communities responded slightly differently to edaphic variables in that they displayed substantial and strong correlations with base saturation, DOC, pH and TP, followed by relatively weak correlations with SOC, TN, clay and CEC. Variables that strongly correlated with the overall microbial community composition were also strongly associated with those of certain bacterial/archaeal phyla and fungal classes as revealed by Spearman R values (Table S4). Among the most abundant bacterial/archaeal phyla, Acidobacteria were negatively correlated with soil pH, base saturation, CEC and TN and were positively correlated with DOC. Actinobacteria and Deltaproteobacteria were positively associated with soil pH, base saturation, CEC and TN.

Fig. 2

Nonmetric multidimensional scaling (NMDS) based on Bray-Curtis dissimilarity displaying sites and significant environmental vectors for bacterial/archaeal (A) and fungal (B) community composition in six soils from three land uses (C: cropland; G: grassland; F: forest) on two bedrock forms (S: silicate; L: limestone). Circles indicate silicate soils, and triangles limestone soils. Black represents cropland soils, green grassland soils, and red forest soils. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 2

Spearman rank correlations (R values) of single edaphic variables with matrices of microbial community composition, extracellular enzyme activities and soil processes based on Mantel tests.

Edaphic variables	Bacterial/archaeal composition	fungal composition	Extracellular enzyme patterns[a]	Soil processes[b]
pH	0.74**	0.44**	0.54**	0.38*
Base saturation	0.78**	0.53**	0.50**	0.39*
CEC	0.51**	0.17*	0.12	0.33**
Clay	0.30**	0.25*	0.10	0.01
SOC	0.32*	0.36*	0.10	0.31*
TN	0.44**	0.27*	0.20	0.21
TP	0.24*	0.43**	0.20	−0.08
DOC	0.57**	0.51**	0.36*	0.21

qGlucosidase, qPhenoloxidase, qAminopeptidase, qPhosphatase.

CUE, qGrowth, qCO2, NUE, qGross protein depolymerization, qGross mineralization, qGross nitrification. Significance levels: ***: P < 0.001; **: P < 0.01; *: P < 0.05.

Relationships between single or combined edaphic variables with microbial communities, soil processes and extracellular enzyme patterns

Based on Mantel tests, soil pH, base saturation and DOC showed significant correlations with the extracellular enzyme matrix while TN, TP, CEC, SOC and clay content were not correlated with enzyme patterns (Table 2). With respect to the soil process matrix, we found substantial correlations with base saturation, pH, CEC and SOC, but no correlations with TN, DOC, clay and TP.

When we combined the significant edaphic variables as shown in Table 2 into an edaphic matrix, the matrix was significantly correlated with both bacterial/archaeal and fungal communities based on Mantel tests (Fig. 3). A stronger effect of the edaphic matrix was observed on bacterial/archaeal communities than on fungal communities. Moreover, according to partial Mantel tests, the edaphic matrix showed a direct correlation with the soil C process matrix and the matrix of all measured soil processes but not with the soil N process matrix or enzyme patterns. The bacterial/archaeal community was significantly correlated with soil N process matrix, which consequently led to a significant association with the combined matrix of all soil processes. Although the fungal community was substantially associated with NUE and qPhosphatase activities (Table S5), it was not correlated with either matrix of soil C or N processes or of enzyme patterns. Neither bacterial/archaeal nor fungal communities were significantly correlated with enzyme patterns. No correlations between enzyme patterns and soil C or N processes were observed.

Fig. 3

Spearman rank correlations (R values) between matrices of edaphic variables, microbial community composition, soil C and N processes and soil enzyme activities based on Mantel tests (red) and partial Mantel tests (black). EDV: edaphic variables including soil pH, clay content, base saturation, CEC, SOC, TN, TP and DOC. BCC: bacterial/archaeal OTU composition. FCC: fungal OTU composition. EEA: extracellular enzyme activities normalized to MBC including qGlucosidase, qPhenoloxidase, qAminopeptidase, and qPhosphatase. C-Process: microbial CUE, qGrowth and qCO2. N-Process: microbial NUE, qGross protein depolymerization, qGross mineralization and qGross nitrification. Soil-Process: C-Process + N-Process. Grey values show non-significant correlations. Significance levels: ***: P < 0.001; **: P < 0.01; *: P < 0.05. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Links between microbial diversity and edaphic factors, soil enzymes and soil processes

Based on one-way ANOVA results, bacterial/archaeal and fungal richness were not affected by land use type, but was significantly lower in silicate soils compared to limestone soils (Chao1 index, P < 0.05). Both bacterial/archaeal (P < 0.01) and fungal (P < 0.001) Shannon diversity were significantly lower in forest soils compared to cropland and grassland soils, but were not influenced by bedrock forms. The richness and diversity of bacterial/archaeal and fungal communities showed substantial correlations with the edaphic matrix based on Mantel tests (Table S6A). Moreover, the richness and diversity of bacterial/archaeal and fungal communities exhibited strong positive associations with soil pH and CEC, and fungal diversity was negatively correlated with DOC content (Table S6B). The richness and diversity of bacterial/archaeal and fungal communities also showed strong negative relationships with extracellular enzyme activities. The richness and diversity of bacterial/archaeal and fungal communities exhibited no correlations with the soil C process matrix but strong correlations with the soil N process matrix. Strong negative associations between bacterial/archaeal diversity and qGross protein depolymerization rate, and between fungal richness and qGrowth and qGross protein depolymerization rate were found according to Spearman rank correlations.

Discussion

Edaphic factors shape microbial community composition and diversity in soils differing in land use and bedrock

The heterogeneous nature of the soil (micro) environment is thought to maintain highly diverse microbial communities (Fierer, 2017). Here, we explored the relationship between a range of edaphic variables and microbial community composition and diversity across three different land uses and two bedrocks. When considering single edaphic factors, most physicochemical properties (e.g., soil reaction, nutrients and texture; but not SOC and TP) showed a stronger effect on bacterial/archaeal community composition than on that of fungi. One possible explanation might be that bacteria/archaea are fostered to better adapted to local edaphic conditions than fungi due to their different growth strategies, as fungi may access more soil volume due to their hyphal growth and thereby get access to more substrates and nutrients than bacteria/archaea. However, in native soil environments, edaphic variables co-vary and likely interact to regulate microbial community structure, diversity and function, since soil environments are defined by a combination of edaphic and climatic characteristics that microorganisms must adapt to in synchrony (Schimel and Schaeffer, 2012). Treating the edaphic properties as a matrix allowed us to investigate their combined impact on microbial community composition. Here, with respect to the edaphic matrix, significant correlations with both bacterial/archaeal and fungal communities were observed, with slightly stronger responses of the bacterial/archaeal communities. The combined edaphic effects on both bacterial/archaeal and fungal communities were much stronger than the effects of single edaphic variables, illustrating for the first time on a matrix level that the combination between those edaphic factors strengthened the environmental influence on microbial community composition compared to the effect of individual factors.

Among all edaphic variables, we found strongest Mantel correlations between bacterial/archaeal as well as fungal communities and soil pH and base saturation (Table 2), both of which are strongly positively related to each other according to PCA (Fig. S3). Soil pH has been widely recognized as a key factor influencing microbial community composition (Lauber et al., 2008; Rousk et al., 2010a, 2010b). Among the most abundant bacterial phyla, Acidobacteria and Actinobacteria exhibited a strong, inverse responses to soil pH (Table S4), corroborating the fact that members of Acidobacteria (e.g. Subdivision 1 and 3) tend to become more prominent at mildly acidic pH (Eichorst et al., 2007; Foesel et al., 2014; Jones et al., 2009; Sait et al., 2006). Actinobacteria were reported to thrive in soils with neutral pH and to grow best between pH 6 to 9 (Barka et al., 2016), as supported by the observed strong positive correlation between Actinobacteria and soil pH in this study. Soil pH was found to exert different or even contrasting effects on bacterial and fungal communities, i.e. low soil pH was found to decrease bacterial growth while to increase fungal growth (Rousk et al., 2009), which could potentially alter the microbial community structure by favoring low-pH adapted or acidophilic microorganisms. Here, we found that both bacterial/archaeal and fungal communities were affected by soil pH, but the bacterial/archaeal community was more strongly influenced by pH than that of fungi, which might be due to relatively narrow optimal pH ranges for bacterial growth but wide pH ranges for fungal growth (Rousk et al., 2010a). Despite the direct influence of soil pH on microbial community structure, soil pH can also shape microbial communities indirectly by other co-varying factors such as nutrient availability and organic C content (Rousk et al., 2010a). Additionally, base saturation, representing the percent of the cation-exchange sites occupied by basic cations such as Ca2+, Mg2+, Na+, and K+, was significantly correlated with both bacterial/archaeal and fungal communities, indicating that base saturation is a global variable that co-explains microbial community dissimilarity. Cations such as manganese (Mn) was found to shape microbial community composition independent of pH (Whalen et al., 2018). The significant impact of base saturation on microbial community composition, however, is most likely explained by the co-variance of soil pH and base saturation in this study, indicating that soils with higher base saturation typically have higher soil pH and generally are more fertile. Moreover, low soil pH might lead to the higher solubility of SOM and altered composition of dissolved organic matter in soils (Curtin et al., 2016), which could trigger changes in energy and nutrient availability for microorganisms and thereby affect microbial abundance and composition. Consistent with previous work (Li et al., 2012; Xu et al., 2015), we also found a significant correlation between DOC and bacterial/archaeal community composition, which was mainly explained by the positive correlation between DOC and Acidobacteria. The reported negative responses of Acidobacteria to increased available organic C suggested that members of this phylum are oligotrophic bacteria (Fierer et al., 2007); however, our results suggest that not necessarily all Acidobacteria are oligotrophs, corroborating the fact that some of the Acidobacteria isolates could grow in higher C concentrations (Kielak et al., 2016; Navarrete et al., 2015).

As indispensable energy and nutrient source for microorganisms, soil organic matter content (as represented by SOC or by soil TN) was reported to play an important role in shaping microbial communities (Burns et al., 2016; Drenovsky et al., 2004). For instance, organic C and N amendment experiments revealed significant changes in microbial PLFA composition and in fungal: bacterial ratios (Drenovsky et al., 2004; Ng et al., 2014; Zhou et al., 2017). The relative abundances of Actinobacteria were reported to increase with soil C and N pool size (Li et al., 2014), in accordance with our finding that Actinobacteria showed positive correlations with SOC and TN. Another dominant microbial phylum, Firmicutes (copiotrophic), was found to be negatively associated with SOC. This contradicts the typically observed positive relationship of Firmicutes with soil C content (Ling et al., 2017; Tsiknia et al., 2014), which may be due to differences in the quality and accessibility of SOC (not specifically measured in this study). If there was less biodegradable or bioaccessible SOC, more SOC will not necessarily lead to a greater abundance of copiotrophic microbial communities including Firmicutes. The Deltaproteobacteria and Acidobacteria were also strongly correlated with TN, corroborating with previous research (Ling et al., 2017; Zhang et al., 2013). There is no general agreement of the effects of phosphorus (P) - another crucial nutrient - on microbial community composition, as negative, neutral or positive effects of P were found on soil microbes in terrestrial ecosystems (DeForest et al., 2012; Huang et al., 2016; Liu et al., 2012). In this study, TP was found to be the most crucial edaphic factor in explaining dissimilarities in fungal communities (NMDS), in line with the reported important role of P in structuring soil fungal communities in P addition experiments (He et al., 2016). TP is typically less in forests than in managed ecosystems due to fertilization. We found higher abundances of Eurotiomycetes and Agaricomycetes in forest soils as compared to cropland and grassland soils (both P < 0.05), which was due to negative correlations between TP and the relative abundance of these two fungal classes. Based on SIMPER analysis (Table S3B), Eurotiomycetes and Agaricomycetes accounted for 48.1% of the overall dissimilarity between forest soils and cropland soils, and explained 52.2% dissimilarity between forest soils and grassland soils. Therefore, TP might be an important driver in structuring soil fungal communities across land uses, though the exact mechanism currently remains elusive.

Despite the influence of land use and bedrock on microbial community composition, we also observed strong effects of bedrock on microbial richness and effects of land use on microbial diversity, which were likely due to the influence of bedrock and land use on soil pH and base saturation (Table 1). Soil pH has a strong impact on microbial diversity across different spatial scales and soil types (Lauber et al., 2009; Rousk et al., 2010b; Zhalnina et al., 2015), corroborating our results that pH was positively correlated with both bacterial/archaeal and fungal diversity (Table S6B). The strong positive correlation between base saturation and microbial diversity is likely due to the strong association between base saturation and soil pH (Fig. S3). The impact of C, N and P on microbial diversity was not consistent in previous observational and experimental studies (Leff et al., 2015; Li et al., 2012; Ling et al., 2017; Lopez-Fernandez et al., 2018). Microbial diversity was negatively associated with DOC content here, in line with the previous results obtained at laboratory and field scales (Li et al., 2012), although positive correlations between DOC and microbial diversity were also reported previously (Lopez-Fernandez et al., 2018). A possible explanation might be grounded in the effect of DOC on microbial community composition. For example, the observed positive correlations between DOC content and the dominant bacterial phylum Acidobacteria and fungal class Eurotiomycetes (Table S4) might lead to less competitiveness and less influence by other bacteria, archaea and fungi in the studied soils, and thus may result in negative correlations between microbial diversity and DOC content. There were no significant relationships between microbial diversity and soil N and P content in this study, which were different from previous research that showed negative correlations between microbial diversity and N and P content (Leff et al., 2015; Ling et al., 2017). Inconsistent with previous studies (Lynn et al., 2017; Ma et al., 2016), we did not find a significant association between microbial diversity and clay content and CEC. Therefore the significant combined effect of the edaphic matrix on microbial diversity demonstrated by Mantel tests (Table S6A) is likely induced by the effect of individual edaphic parameters including soil pH, base saturation and DOC content on microbial diversity.

The influence of edaphic variables, microbial community composition and diversity on soil C and N processes

The soil environmental variables are generally regarded as good predictors of soil C and N process rates (Graham et al., 2016). For example, soil pH is often positively linked with substrate and nutrient availability (Mccauley et al., 2017) and is expected to affect soil microbial C and N processes. In this study, we found that soil microbial C and N processes were significantly correlated with individual edaphic variables (Table S7). For instance, we observed a strong negative correlation between soil pH and qGrowth, likely due to the negative influence of pH on DOC concentrations (Pearson R = −0.73, P < 0.001). DOC represents a major labile C and energy source for microbes and was found to positively affect qGrowth in soils (Zheng et al., 2019). It is likely that the strong effects of soil pH, base saturation, CEC and DOC on qGrowth ultimately led to a substantial correlation between the edaphic matrix and the soil C process matrix (Fig. 3). Single edaphic variables also exhibited strong connections to some of the soil N processes or extracellular enzyme activities, e.g. SOC and TN content showed negative connections with gross protein depolymerization and gross mineralization, and TP, pH and base saturation were all negatively associated with qPhenoloxidase and qPhosphatase activities. However, when considered as a matrix, edaphic properties showed no significant effect on the soil N process matrix or extracellular enzyme patterns, indicating no or very weak influences of the combined edaphic properties on soil N processes and soil enzyme patterns. When incorporating all measured soil C and N processes into one soil process matrix (Fig. 3), the connection of the edaphic matrix to this merged soil process matrix became stronger (R = 0.36, P < 0.05) than the connection of the edaphic matrix to the soil C process matrix (R = 0.30, P < 0.05) or to the soil N process matrix (R = 0.23, P > 0.05) individually. This demonstrates that the influence of edaphic properties on soil processes strengthened when more processes were incorporated into the soil process matrix, i.e. the more multifunctional the consideration of soil processes became. This again illustrates the importance of investigating the effects of multiple edaphic factors on multiple soil functions instead of only studying the relation between single soil parameters and single soil processes.

Microbial community composition has been variably demonstrated to affect microbial processes (Becker et al., 2017; de Menezes et al., 2017; Graham et al., 2016). Here partial Mantel tests showed that bacterial/archaeal community composition was significantly affecting single soil C or N processes, i.e. microbial growth, microbial NUE, gross protein depolymerization and gross N mineralization rates (Table S5), corroborating findings of previous studies on soil respiration, net N mineralization and denitrification (Colman and Schimel, 2013; Li et al., 2015; Zhou et al., 2011). In terms of the dominant archaeal phylum Thaumarchaeota, members of which were found to oxidize ammonia aerobically and contribute to the soil nitrification process (Brochier-Armanet et al., 2012; Pester et al., 2011), showed no correlation with nitrification rates in this study. This unambiguously demonstrates the influences of bacterial/archaeal community composition on specific soil processes. Although the bacterial/archaeal community composition was not significantly correlated with the soil C process matrix, its significant correlation with qGrowth highlights that some carbon transformation processes are inherently linked to bacterial/archaeal community composition. Here microbial respiration (as represented by qCO2) was not correlated with microbial community composition (Table S5), in accordance with other short-term studies (Barnard et al., 2015; Leff et al., 2012; Placella et al., 2012). Only few studies recorded a relationship between microbial respiration and specific bacterial lineages (Che et al., 2016; Fierer et al., 2007; Orr et al., 2015), but the results were not consistent across studies. Here, no significant correlations between microbial respiration and specific bacterial/archaeal or fungal lineages were detected (Table S8). Almost no studies tested the relationships between microbial respiration and archaeal or fungal lineages (Che et al., 2016), and we did not find substantial correlations between microbial respiration and archaeal or fungal lineages in this study. Additionally, bacterial/archaeal community composition exhibited a strong regulatory effect on soil N processes (matrix-level) likely due to its significant correlations with single soil N processes (Table S5). Similarly, fungal community composition showed no direct influence on combined C or N process matrices possibly due to its weak influences at the single process levels. The reason for this was possibly the low activity of the fungal community, which was not specifically measured in this study, compared to the bacterial/archaeal community across the studied soils. The observed significant correlation between gross N mineralization and bacterial community composition was likely due to the strong correlations between gross N mineralization rate and two dominant bacterial lineages, i.e. Acidobacteria and Deltaproteo-bacteria. Likewise, no significant association between fungal community composition and gross N mineralization rate was found since most abundant fungal classes were not correlated with this process here. Although not all of the single processes were associated with bacterial/archaeal community composition, we still found a strong control of the bacterial/archaeal community composition on soil C and N processes, almost rivaling the direct edaphic effects on merged soil processes. This again highlights that edaphic and microbial controls on soil processes strengthen when considering soil processes in a multifunctional context and edaphic variables not in an isolated but combined form.

Despite the significance of microbial community composition in regulating multiple soil processes, soil microbial diversity also plays a pivotal role in maintaining ecosystem multifunctionality (Delgado-Baquerizo et al., 2017b, 2016; Wagg et al., 2014), corroborating the strong correlations between bacterial/archaeal/fungal richness/diversity and the soil process matrix observed here (Table S6A). Interestingly, no effects of microbial richness or diversity on the soil C process matrix or on microbial respiration or CUE were detected, in accordance with previous field study that showed no relation between basal respiration with bacterial richness (Delgado-Baquerizo et al., 2017b). However, a strong correlation between fungal richness and qGrowth was observed, illustrating that microbial richness may play a role in influencing microbial growth. In contrast, we found strong associations between microbial richness/diversity and the soil N process matrix (Table S6A), likely induced by strong effects of microbial richness/diversity on gross protein depolymerization (Table S6B), which has been rarely reported previously. In general, compared to bacterial/archaeal richness and diversity, fungal richness and diversity seemed to play stronger roles in regulating soil N processes as demonstrated by partial Mantel results.

Except for direct influences, microbial communities could also affect soil processes by regulating extracellular enzyme levels. Extracellular enzymes are essential for organic matter decomposition and are primarily produced by fungi and bacteria in soils. Fungi are generally considered to process recalcitrant C and N-poor substrates while bacteria are thought to be more responsive to labile substrates (Treseder et al., 2016; Xu et al., 2015). Moreover, fungi are thought to possess a greater capacity to produce extracellular enzymes for decomposition of complex plant organic matter than bacteria, and intermediate decomposition products by fungi can provide labile resources for bacteria (Romaní et al., 2006). In this study, soil enzyme patterns were not significantly correlated with either the bacterial/archaeal or the fungal communities, in contrast to previous studies showing that microbial community composition shaped enzyme patterns (Gallo et al., 2004; Schnecker et al., 2015; Waldrop et al., 2000). A possible explanation for our finding is that here we only included four extracellular enzymes involved in the decomposition of soil organic C, N and P into the matrix. It is therefore possible that having data on more divergent extracellular enzymes may lead to a stronger impact of bacteria/archaea and fungi on enzyme patterns. Nevertheless, we still found significant correlations of qPhenoloxidase and qPhosphatase with bacterial/archaeal community composition, as well as a substantial effect of fungal community composition on qPhosphatase, in accordance with the reported influence of microbial community composition on soil enzyme activities (Talbot et al., 2013; Waldrop and Firestone, 2006). We also observed substantial associations between microbial richness/diversity and soil enzyme patterns and single soil enzyme activities (Tables S6A and B), in line with findings from previous soil multifunctionality studies (Delgado-Baquerizo et al., 2017b, 2017a). Additionally, soil enzyme activities are influenced by environmental factors such as soil pH, soil texture, and substrate and nutrient availability (Acosta-Martínez et al., 2007; Turner, 2010). We also found strong correlations between edaphic factors (e.g. soil pH, base saturation, TP, DOC) and some of the extracellular enzyme activities as well as soil enzyme patterns; however, no significant correlation was found between the edaphic matrix and soil enzyme patterns. This finally indicates that the interactions between edaphic variables could possibly offset their overall effects on enzyme activities compared to the effects of individual edaphic variables.

Conclusions

We found that soil bacterial/archaeal and fungal communities were shaped by similar edaphic variables though slightly differing in strength. Not only single edaphic variables but also the combined multivariate representation of the soil environment strongly affected microbial community composition and diversity. Strong relations between bacterial/archaeal community composition and soil processes were also found, with lesser effects of fungal community composition likely due to low activity. Both bacterial/archaeal and fungal diversity showed strong relations with soil N processes but not with soil C processes. Soil enzyme patterns were not affected by the multivariate soil environment or by microbial community composition, but were shaped by microbial richness and diversity. Moreover, stronger effects of the soil environment on combined soil processes than on soil C processes or N processes individually when studied in multidimensional and multi-functional space were apparent. The limitation of our study is that the results of observational studies are correlative and potentially non-causative, but still we provide useful information on how microbial community composition and diversity relate to the soil environment and to soil multifunctionality under “real world” conditions. Combined manipulative and observational approaches are suggested in future studies of environment-microbial community structure-function interactions.

In conclusion, this study adds to an integrated understanding of using microbial community structure (i.e. composition, richness and diversity) to improve predictions of multiple soil functions (i.e. enzyme patterns, C and N cycling), and confirms the significance of treating the soil environment as an integrity to predict soil multifunctionality.

Supplementary Material

Supplementary data to this article can be found online at https://doi.org/10.1016/j.soilbio.2019.107521.