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).