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