2. METHODS
2.1 Study site
This study was conducted in the Mu Us Desert, south of the Ordos Plateau
of China, with an area of approximately 4 × 104km2 (N38°55’, E 109°53’). The average altitude is
between 1100 and 1300 m. The climate is semi-arid (rainfall: 440
mm/year), and the mean annual temperature is 7.7°C. The soil is mainly
composed of sand and loessal soil. The vegetation is mainly composed of
xerophytes such as Artemisia desertorum , Hedysarum laeve ,Psammochloa villosa , and Hedysarum scoparium .
Based on standardized vegetation coverage criteria (Guo-dong, 2004),
sites at three different desertification stages were selected as
sampling sites using a Global Positioning System (GPS) unit, including
potential desertification (PD), moderate desertification (MD), and
severe desertification (SD) (Table S1). In July 2018, we built ten
replicate plots (20 m × 20 m) for PD and MD due to the high variability
of vegetation cover and edatope, respectively. We built five replicate
plots (20 m × 20 m) for the SD stage due to the low and indistinctive
vegetation cover and edatope. These 25
plots were randomly selected at each
site, with an interval of at least 200 m. Nine subsamples were obtained
by S -type sampling and mixed into one sample from each soil layer
(0 - 10, 10 - 20, 20 - 30, 30 - 60, and 60 - 100 cm) using a soil auger
with a diameter of 5 cm. Each soil sample (125 samples in total) was
fully mixed and divided into two parts. One part was stored at -80°C
until DNA extraction, and the other part was air dried for
physicochemical analysis.
2.2 Edaphic variables and functional
genes
Soil physicochemical analysis of organic carbon (OC),
total nitrogen (TN), nitrate
(NO3−-N), ammonium
(NH4+-N), total phosphorus (TP),
available phosphorus (AP), pH, electrical conductivity (EC), bulk
density (BD) and soil moisture was
conducted by using standard methods described elsewhere (Bao, 2000;
Lozano et al., 2014; Wang et al., 2017). Microbial DNA was extracted
from 0.6~1.0 g of each desert soil (125 samples) by
using the MP Fast DNA Spin Kit for Soil (MP Biomedicals, Santa Ana, CA,
USA) according to the manufacturer’s protocol. Total
Ca2+ was measured by flame atomic absorption
spectrophotometry (AAS), and Zn2+ was measured
sequentially by graphite furnace atomic AAS on a Hitachi Z‐2000
(Hitachi, Tokyo, Japan).
The abundance of functional genes associated with key microbial
processes was quantified to assess multifunctionality. The abundances of
functional genes associated with the C fixation process (CBBLgene), N fixation process (nifH gene), P dissolution process
(phoD gene), N mineralization process (apr and chiAgenes ), NH4+-N transformation
process (AOA and AOB genes),
NH4+-N loss process (anammoxgene), and NO3−-N loss process
(denitrification, including nirK , nirS , qnorB , andnosZ genes) were quantified by
a real-time PCR system (QuantStudio
6 Flex, Thermo Fisher, Singapore City, Singapore). These genes were
described elsewhere (Levy-Booth et al., 2014; Li et al., 2018; Wang et
al., 2017). Detailed information on the gene primers is given in Table
S2. Reaction mixtures (20 μL) contained 10 μL of SYBR® Premix ExTaq™
(Dining), 0.5 μL of DNA, 0.5 μL of primers (forward and reverse) (20
μM), and 8.5 μL of sterile distilled water. At the same time, we tested
the amplification specificity by constructing a dissolution curve.
Amplicon libraries were constructed according to the dual indexing
strategy. The V4-V5 hypervariable regions were targeted using primer
pair 515F (GTGCCAGCMGCCGCGGTAA)/907R (CCGTCAATTCCTTTGAGTTT) for the
bacterial 16S rRNA gene and primer pair Arch519F (CAGCCGCCGCGGT
AA)/Arch915R (GTGCTCCCCCGCCAATTCCT) for the archaeal 6S rRNA gene. The
fungal ITS1 gene was amplified by the primer pair ITS5-1737F
(GGAAGTAAAAGTCGTAACAAGG)/ITS2-2043R (GCTGCGTTCTTCATCGATGC). All the
samples were sequenced on the Illumina MiSeq (300-bp paired-end reads)
platform (Illumina Inc., San Diego, USA) at Personal Biotechnology Co.,
Ltd. (Shanghai, China). As a standard data analysis
method, the acquired sequences were filtered for quality, and any
chimeric sequences were removed with the Quantitative Insights into
Microbial Ecology (QIIME, v 1.8.0, http://qiime.org/) tool. These
original sequences were filtered for quality and assigned to operational
taxonomic units (OTUs) (DeSantis et al., 2006). The OTUs with fewer than
two sequences were first removed, and then the OTUs were classified
within the SILVA database (release 132, http://www.arb-silva.de) for
bacteria and archaea and the UNITE database (release 7.2,
http://unite.ut.ee/index.php) for fungi. The soil microbiome dataset has
been deposited in the NCBI Sequence Read Archive under accession number
PRJNA541211.
2.3 Assessing
multifunctionality
Multifunctionality has been broadly defined as the simultaneous
provisioning of multiple functions and services (Manning et al., 2018).
These processes involve, among others, soil nutrient turnover (i.e., OC
and TN contents, N and P availability, and N mineralization and loss)
and primary production (i.e., net primary productivity) (Manning et al.,
2018; Wagg et al., 2019). Compared to superficial soils, deep soils are
mainly anaerobic and low-temperature environments. Although maximal
nutrient turnover rates can be measured using indoor culture
experiments, it may be largely inappropriate to adopt indoor culture
experiments to measure nutrient turnover rates in deep soils. The
abundance of functional genes has been identified as the predominant
variable to predict potential nutrient turnover rates (Attard et al.,
2011; Graham et al., 2014; Petersen et al., 2012; Wang et al., 2017).
Thus, to summarize the changes in the overall functioning of deep soils,
the abundance of functional genes affiliated with the nutrient turnover
process was used to indirectly characterize their potential functioning
capacity (de Vries et al., 2018; Graham et al., 2014; Li et al., 2018;
Zhi et al., 2015). In desert ecosystems, the desertification process is
characterized by the long-term deterioration of soil multifunctionality
and appearance of environmental stress gradients for soil microorganisms
(Fierer, 2017; Ravi et al., 2010). Thus, functional gene abundance is
likely to be more accurate than other metrics as a long-term proxy and
predictor of the in-site functioning capacity in deep soils (Li et al.,
2018; Wang et al., 2019b). Based on well-known methods, the abundance of
functional genes associated with functional processes (i.e., C, N and P
transformation) was organized into a multifunctionality index to obtain
and summarize a multidimensional multifunctionality index in deep soils.
In total, we quantified multifunctionality using related variables as
follows: C functions (CBBL ), N functions (TN,
NH4+-N,
NO3--N, nifH , apr ,chiA, AOA , AOB, anammox , nirK ,nirS , qnorB , and nosZ ), P functions (TP, AP, andphoD ), Ca2+ and Zn2+.
There are several approaches to quantify multifunctionality (i.e.,
single functions, averaging, and the threshold approach) (Manning et
al., 2018). Because each approach possesses different strengths and
weaknesses, we adopted two distinct approaches (single functions and the
averaging approach) to quantify and assess multifunctionality. We
normalized and standardized each nutrient cycling variable and
functional gene abundance using Z-score transformation (Manning et al.,
2018). These standardized ecosystem functions were then averaged to
obtain an overall multifunctionality index and single-function index
(i.e., C process, N process, and P process). This index is broadly
adopted in multifunctionality research and provides a straightforward
measure of the ability of different communities to sustain
multifunctionality (Manning et al., 2018; Wagg et al., 2019). In
addition, functional gene abundances for
NH4+-N and
NO3−-N loss were inverted, as higher
abundances of these genes are undesired functions (i.e., a lower
abundance of N loss indicates high N retention and better ecosystem
functioning, while a higher abundance indicates dysfunction and higher N
loss).
2.4 Data analysis
Z-score values of ecosystem functions and related variables were
calculated using SPSS (IBM SPSS Statistics 25.0). The most abundant
phyla and alpha-diversity indexes (sobs, ACE, Shannon, and phylogenetic
diversity) of soil bacteria, archaea, and fungi were calculated in the R
environment (v3.5.3; http://www.r-project.org/) using MOTHUR R package
(https://mothur.org/wiki/calculators/), and their variations with
increasing soil depth were determined using linear least-squares
regression models (OriginLab Corporation, One Roundhouse Plaza, Suite
303, Northampton, MA 01060, United States). The beta-diversity of the
bacteria, archaea, and fungi was quantified by using two axes from a
nonmetric multidimensional scaling (NMDS) analysis of Bray-Curtis
dissimilarities within the OTU community matrix (Daum & Schenk, 1997).
Similarity values among the sites were examined via the ANOSIM and
ADONIS test (P < 0.001). The vertical spatial decay
relationships were calculated as the linear least-squares regression
relationships between soil depth and microbial community similarity
(Jiao et al., 2018).
Classification random forest (RF) analysis was adopted to identify the
importance of each predictor of soil multifunctionality (Wagg et al.,
2014). In these RF models, soil variables (i.e., depth, pH, EC, and BD)
and the most abundant phyla and alpha-diversity indexes of bacteria,
archaea, and fungi served as predictors of multifunctionality. These
analyses were conducted using the randomForest package in the R
environment (Delgado-Baquerizo et al., 2016). The significance of the
importance of each predictor of multifunctionality was identified by
using the rfPermute package for R
(https://cran.r-project.org/web/packages/rfPermute/index.html).
Structural equation modelling (SEM) was performed to determine the
direct and indirect effects of soil variables (i.e., depth, pH, EC, and
BD) and the most abundant phyla and alpha-diversity of bacteria,
archaea, and fungi on multifunctionality. We used the chi-square
(χ2) test and root MSE of approximation (RMSEA) to
assess SEM fit. The model has a good fit when
χ2/df (degrees of freedom) is low (≤ 2) andP is high (> 0.05) and when RMSEA is low (≤ 0.05)
and P is high (> 0.05) (Schermelleh-Engel et al.,
2003). Furthermore, we evaluated the fit of the models using the
Bollen-Stine bootstrap test. Each SEM was fit using AMOS 20.0
(International Business Machines Corporation, Armonk, NY, United
States).