2 MATERIALS AND METHODS
2.1 Study area and soil sampling
The sampling site was located in Puding county, Anshun city, southwest Guizhou Province, China (26°22′–26°26′N and 105°34′–105°41′E). The region is characterized by a subtropical monsoon climate, with annual average sunshine of 1165 h, annual average temperature of 15.1°C, frost-free period of 301 d, and annual average rainfall of 1378 mm. The karst landforms in the territory are widely developed and complete. The soil type is limestone soil according to Chinese soil general classification, which is similar to Mollic Inceptisols according to USDA Soil Taxonomy (Soil Survey Staff, 2010) .
Four ecosystems covering different land-use intensity levels were selected according to the degree of anthropogenic perturbation: (1) primary forest (PF, undisturbed natural vegetation, dominant plant species: Lithocarpus confinis, Itea yunnanensis and Carpinus pubescens); (2) secondary forest (SF, without cultivation and fertilization >15 y, dominant plant species: Rhus chinensis, Populus adenopoda Maxxim and Rosa cymosa); (3) abandoned land (AL, without cultivation and fertilization <6 y, dominant plant species: Conyza canadensis, Imperata cylindrica and Artemisia dubia); (4) cultivated land (CL, cultivation and fertilization with urea/compound fertilizer with some excreta, dominant plant species: Brassica napus and Zea mays). In April 2017, soil samples (n = 6) from 0-10 cm depth were collected using a trowel in each land use type. In brief, Typical plots (10 m×20 m) of each land-use type were sampled from six established sites at six slope positions at least 30 m apart along one altitudinal transect. The 4 replicate soil samples taken at each site were composited to provide one representative sample for each site. The fresh soil samples were sieved to 2 mm and visible roots and stones were removed. For each samples, one portion was air-dried and used for soil chemistry and mineral content analyses, while the other portion was immediately frozen at -80 °C for subsequent microbial analyses.
2.2 Soil physiochemical analyses
Soil pH was determined in 1:2.5 (weight:volume) suspensions by ULTRAMETERII™6PFCE (Myron L, CA, USA). The samples used for elemental analysis of iron (Fe), calcium (Ca), magnesium (Mg), aluminium (Al) and sulphur (S) were prepared by dissolution in HCl-HNO3-HF-HClO4 (Zhao et al., 2016) and analyzed using an inductively coupled plasma optical emission spectrometer (ICP-OES, Optima 5300 DV, Perkin Elmer, Waltham, MA, USA).
2.3 Characterizing soil microbial communities
For each sample, DNA was extracted from 0.25 g of frozen soil (Power soil DNA Isolation Kit, MoBio Laboratories, Carlsbad, CA, USA), and quantified by spectrophotometry (Nanodrop Technologies, USA). The V3-V4 hypervariable regions of bacterial 16S ribosomal RNA (rRNA) encoding gene were amplified using the primers 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′- GGACTACHVGGGTWTCTAAT-3′). PCR was performed in triplicate using the following thermal program: 95 °C for 3 min, followed by 27 cycles at 95 °C for 30 s, 55 °C for 30 s, 72 °C for 45 s and a final extension at 72 °C for 10 min. The PCR products were subjected to electrophoresis on 2% agarose gels for detection, purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA), and then quantified using QuantiFluor™-ST (Promega, USA).
Illumina next-generation amplicon sequencing was used to sequence the samples (Illumina MiSeq platform, Illumina, San Diego, CA, USA) . Raw sequences >200 bp with an average quality score >20 and without ambiguous base calls were quality processed, using the Quantitative Insights into Microbial Ecology (qiime) pipeline (version 1.17). Operational taxonomic units (OTUs) were picked at 97% sequence similarity. The resultant OTU abundance tables from these analyses were rarefied to an even number of sequences per samples to ensure equal sampling depth (14371 for 16S rDNA)
2.4 Assessing multiple ecosystem functions and the resistance of multifunctionality to human disturbance
On all soil samples, thirteen functions related to C, N and P cycling were measured: activity of extracellular soil enzymes α-1,4-glucosidase (αG), β-1,4-glucosidase (βG), β-xylose (βX), cellulose diglucosidase (CBH), N-acetylglucosaminidase (NAG), leucine aminopeptidase (LAP) and phosphatase (AP), dissolved and total organic carbon (C), total nitrogen (N), extractable ammonium (NH4+) and nitrate (NO3-) concentrations and total P concentration (TP). The total C and N contents, enzyme activities, inorganic nitrogen concentrations and AP and TP contents were obtained using methods as described by Wang et al. (2019) and presented in Table S1. Collectively, these variables constitute a good proxy for the processes contributing to nutrient cycling (Maestre et al., 2012).
We then used the Orwin & Wardle (2004) index (RS) to evaluate the multifunctionality resistance as:
RS=\(1-\frac{2\mid D0\mid}{(C0+\mid D0\mid)}\) (1)
where, D0 is the difference between the natural ecosystem (C0, value of each functional variable of primary forest) and the disturbed ecosystem (P0, secondary forest, abandoned land or cultivated land) soils. This index has the advantage of being: (1) standardized by the control and (2) bound between -1 (lowest resistance) and +1 (maximal resistance) even when extreme values are encountered (Orwin & Wardle 2004). The resistance of each function was calculated independently for each land-use intensity level. In addition, the resistance of the thirteen functions measured was then averaged for each land-use intensity level to obtain a standard index of multifunctionality resistance, which is more integrative than the response of individual functions and better reflects overall ecosystem functioning. The same method has been used to determine indexes of multistability (Durán, et al., 2017) as well as multifunctionality (Delgado-Baquerizo et al., 2016). Finally, the resistance of bacterial taxa also was calculated using equation (1). In this case D0 is the difference between the natural ecosystem (C0, bacterial relative abundance of primary forest) and the disturbed ecosystem (P0, secondary forest, abandoned land or cultivated land) bacterial relative abundance.
2.5 Statistical analyses
Structural equation modelling (SEM; Grace 2006) was used to evaluate the direct and indirect relationships between land-use intensity, soil pH, soil elemental characteristics(Fe, Mg, Ca S and Al) and soil microbial composition (phylum- and class-level predictors ) and ⍺-diversity on multifunctional resistance. The overall goodness-of-fit of the SEM results were evaluated by Chi square test (𝜒2 ; the model has a good fit when 0 ≤ 𝜒2 /df ≤ 2 and 0.05 < P ≤ 1.00), root mean square error of approximation (RMSEA; the model has a good fit when 0 ≤ RMSEA ≤ 0.05 and 0.10 < P ≤ 1.00 –Schermelleh-Engelet al., 2003). Additionally, the standard total effects (STE) was calculated for each variables included in the SEM.
To gain a mechanistic understanding of the drivers of resistance of multifunctionality, within the bacterial community, we conducted a Random Forest classification analysis (Breiman 2001; Delgado-Baquerizo et al., 2016) to identify the best class-level microbial predictors based on the relative abundance of bacterial taxa. We analyzed the importance of soil microbial communities as the predictor of both multifunctional resistance and to individual functions. These analyses were conducted using the random forest package (R Core Team, 2019) of the R statistical software, version 3.0.2. The significance of the importance of each predictor was assessed using the rfPermute package (Archer, 2016).