Warming modulates soil multifunctionality through assembly processes and co-occurrence patterns of arbuscular mycorrhizal fungal communities in drylands

Climate warming poses a threat to the functionality of global dryland ecosystems, while arbuscular mycorrhizal fungal (AMF) communities play a critical role in maintaining ecosystem stability and functioning. However, it remains unclear how the assembly mechanisms and co-occurrence network patterns of AMF communities respond to warming and regulate soil multifunctionality. We conducted a three-year warming experiment (+0.5－1.6°C) using open-top chambers (OTCs) in the Tengger Desert to investigate the assembly processes and co-occurrence network patterns of AMF in the rhizosphere soil of dominant shrubs, Artemisia ordosica and Caragana korshinskii. By measuring 15 soil functional indicators (including nutrient availability, biogeochemical cycles, microbial activity, and microbial productivity), we constructed a soil multifunctionality index and employed structural equation modeling (SEM) to uncover the mechanisms by which AMF communities regulate soil multifunctionality. Our results showed that warming increased the network complexity of AMF communities associated with C. korshinskii by 8%－145%, while decreasing that of the A. ordosica by 14%－80%. Stochastic processes dominated the assembly of AMF communities, with warming reducing the stochasticity of AMF communities in C. korshinskii while increasing it in A. ordosica. Notably, the two species showed distinct functional response pathways: For A. ordosica, warming-induced reduction in soil water content suppressed the mycorrhizal colonization rate (36%－79%) through a "stochastic assembly-network collapse" cascade effect, subsequently leading to decreased soil multifunctionality (54%－172%). In contrast, for C. korshinskii, the soil water reduction caused by warming maintained higher colonization rates (82%－198%) by enhancing AMF network complexity, thereby improving soil multifunctionality (43%－228%). We further proposed a conceptual framework that integrates niche theory into a mechanistic understanding of how stochastic processes and network complexity of AMF communities affect soil multifunctionality under climate warming. Overall, our study suggests that warming modulates soil multifunctionality by affecting the assembly processes and network complexity of AMF communities, exhibiting species-specific responses. This provides crucial theoretical support for understanding the mechanisms underlying functional evolution in drylands and for formulating climate adaptation strategies.

Methods
Study site and field experimental design

Samples were collected at the Shapotou Desert Research and Experiment Station of the Chinese Academy of Sciences (SDRES, CAS; 37°32′ N, 105°02′ E), located on the southeastern edge of the Tengger Desert in North China. The mean annual temperature at an altitude of 1330 m over the past ten years is approximately 8.6°C, ranging from mean average temperatures of -6.9°C in winter (January) to 24.3°C in summer (July). The mean annual precipitation is approximately 186 mm, with 80% falling during the warmer growing season (May－September). Long-term observations from SDRES (1955–2016) show that the average annual temperature has increased by 0.5°C while the average annual precipitation has decreased by 4 mm over the last decade (Li et al., 2021; Li et al., 2023).

The experimental design of this site has been described in detail in previous studies (Li et al., 2021; Li et al., 2023). In brief, the open-top chambers (OTCs) were deployed in the field in 2019, with their lower 50 cm buried underground. Seedlings were transplanted into the OTCs in May 2020, followed by the initiation of official monitoring. These experimental setups were specifically designed to model the temperature increase in the Tengger Desert over the previous decade and the current one. To accomplish this, we systematically adjusted the shape and size of the OTCs to simulate the targeted warming conditions. Ultimately, we selected four OTCs of different sizes from a series of existing devices and used them to simulate four treatments in the experiment. The four different sizes of OTCs are: a large octagonal column of 130 cm × 130 cm × 200 cm (+0.5°C warming); a small octagonal column of 100 cm × 100 cm × 200 cm (+0.9°C warming); a large quadrilateral of 130 cm × 130 cm × 200 cm (+1.2°C warming); a small quadrilateral of 100 cm × 100 cm × 150 cm (+1.6°C warming). In addition, we also set up control groups that were at least 2 meters away from the OTCs to minimize bias. These control plots replicated the OTC subsurface conditions (buried 50 cm deep) and underwent identical pretreatment, installation, and recovery protocols, differing only in the absence of aboveground OTC structures. Air temperature, precipitation, and soil water at 0-5 cm depth in each OTC and the control group were continuously monitored throughout the experiment using an automated weather station (U30 HOBO data logger, Onset Computer Corporation, Bourne, MA, USA). Analysis of the data showed that the mean annual temperature of large, medium, and small OTCs increased by 0.5°C, 0.9°C, 1.2°C, and 1.6°C, respectively, compared to the control.

Soil sampling and property measurement

All field sampling for this study was conducted in mid-September 2023, three years after the implementation of the OTC. The experimental plants, A. ordosica and C. korshinskii, were cultivated as follows: seeds were collected from mature parent plants of the same species in the study region, screened for uniformity in size, shape, and thousand-seed weight, and germinated in a controlled nursery. Seedlings were bare-root transplanted into experimental plots in May 2020 after reaching 7－8 cm in height. Monoculture treatments contained either A. ordosica (1 plant/m²) or C. korshinskii (2 plants/m²), planted in pre-cleared sandy substrate with no pre-existing vegetation. No other plant species were introduced into the OTCs during the experiment.

In each OTC and control treatment, we collected the above-ground parts and below-ground root systems of individually planted A. ordosica and C. korshinskii plants. Four to five plants per treatment were randomly selected for collection, and their root systems were excavated to a depth of 50 cm. Rhizosphere soil collected from the 20－40 cm soil layer of these plants was thoroughly mixed to form a single representative sample. In total, 60 samples (5 treatments × 2 plant species× 6 replicates) were collected. Specifically, during sampling, excess soil was shaken off from the roots by hand, leaving approximately 1 mm of soil still attached to the roots. The roots were then cut into approximately 5 cm segments using sterile tools, sealed in special 50 ml tubes, and immediately placed in liquid nitrogen for quick freezing. These roots were transported to the laboratory to isolate roots and rhizosphere soil using the methods outlined by Zhong et al. (2022). The rhizosphere soil suspension isolated from the roots was stored in a 50 ml Falcon tube and refrigerated at -80°C for subsequent DNA extraction. In addition, the other fresh soil samples were partly present at 4°C to measure microbial biomass and enzyme activity and partly air-dried for soil physicochemical analyses. Roots (diameter < 1 mm) used for analysis of AMF colonization rate were rinsed with distilled water in the laboratory and placed in centrifuge tubes. Both above- and below-ground biomass collected were oven-dried (60°C) to a constant weight.

AMF colonization of roots 

Fine roots (root diameter < 1 mm) were separated from the sieved soil, cut into 1 cm segments, and washed in a 10% (w/v) potassium hydroxide (KOH) solution for 30 min at 90°C (Wu et al., 2023). After rinsing the roots with KOH, the samples were acidified in 2% HCl (v/v) for 30 min and then immersed in 0.05% (w/v) trypan blue for 30 min. Finally, the percentage of root length colonized by AMF, i.e., the mycorrhizal colonization rate (%), was determined using the magnified intersection method (McGonigle et al., 1990).

Amplicon sequencing and bioinformatics analysis 

Following the manufacturer's recommendations, DNA was extracted from 0.25 g of frozen soil for each soil sample using the PowerSoil® DNA Isolation Kit (MoBio Laboratories, Inc., Carlsbad, CA, USA). For Illumina MiSeq sequencing, 18S rRNA gene fragments were amplified by a two-step PCR procedure using GeoA-2 and AML2 primers (Lee et al., 2008) and NS31 and AMDGR primers (Simon et al., 1992; Sato et al., 2005), respectively. The first PCR reaction used a 25 µl reaction mixture with 2.5 µl DNA template, 1.25 µl GeoA-2 and AML2 primers (10 µM each), 25 µl 2 × SYBR Premix Ex Taq (Takara, Japan), and 20 µl sterile water. The first-step PCR product was then diluted 10-fold to serve as the template DNA for the second PCR reaction. The primer pairs in the second PCR were supplemented with the 454 pyrosequencing adapters and the 8-bp barcode for multiplexing. PCR reaction conditions were as follows: 95℃ for 15 min; five cycles at 42℃ for 30 s, 72℃ for 90 s, 92℃ for 45 s; 35 (first reaction) / 20 (second reaction) cycles at 65℃ for 30 s, 72℃ for 90 s, 92℃ for 45 s; 65℃ for 30 s and with a final extension at 72℃ for 10 mins. PCR products were separated on a 2% agarose gel by electrophoresis, purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA), and quantified using the Quantus™ Fluorometer (Promega, USA). The purified PCR products were pooled at equimolar concentrations and subjected to paired-end sequencing on the Illumina MiSeq PE300 sequencing platform (Illumina, San Diego, USA) by standard protocols, as conducted by Majorbio Bio-Pharm Technology Co., Ltd (Shanghai, China).

Raw sequences were filtered using Quantitative Insights into Microbial Ecology (QIIME) 1.9.0-dev pipeline (http://www.qiime.org) (Caporaso et al., 2010). Briefly, sequences with a quality score below 30, less than 200 bp in length, containing ambiguous bases, and no valid primer or barcode were excluded from further analysis. The presence of probable chimeras was examined using UCHIME by the MaarjAM Glomeromycotina database as a reference, and all possible chimeric sequences were removed from the dataset. Non-chimeric sequences were clustered into different operational taxonomic units (OTUs) at a 97% identity threshold using USEARCH v.8.0 based on the UPARSE pipeline (Edgar, 2013). Representative sequences from each OTU were selected and searched against the National Center for Biotechnology Information (NCBI) nucleotide database and the MaarjAM 18S rRNA gene database. The number of sequences per sample was normalized to the minimum sample size using the SUB. SAMPLE command in MOTHUR to account for the effect of different sequence numbers on the analysis of AMF communities. All raw sequencing data have been deposited in the NCBI Sequence Read Archive (SRA) under study accession number PRJNA988244.

Assessment of soil multifunctionality 

We evaluated soil multifunctionality using 15 variables related to nutrient availability, biogeochemical cycles, microbial activity, and microbial productivity (Hu et al., 2021; Qiu et al., 2021). Specifically, these variables included total organic carbon (TOC), total nitrogen (TN), total phosphorus (TP), nitrate nitrogen (NO3- - N), ammonium nitrogen (NH4+ - N), available nitrogen (AN), dissolved organic carbon (DOC), dissolved organic nitrogen (DON), β-1,4-glucosidase (βG), 1,4-β-Dcellobiohydrolase (CBH), β-1,4-N-acetylglucosaminidase (NAG), L-leucine aminopeptidase (LAP), alkaline phosphatase (ALP), and microbial biomass carbon (MBC), and nitrogen (MBN) contents. In addition to obtaining quantitative multifunctional indices for each warming treatment, we first standardized all individual soil functions using z-score transformations (Hu et al., 2021). These standardized z-scores were then averaged to determine the soil multifunctional indices under each warming treatment (Maestre et al., 2012). Extracellular soil enzyme activities such as β-1,4-glucosidase, cellobiose hydrolysis enzyme, β-1,4-N-acetylglucosaminidase, L-leucine aminopeptidase, and alkaline phosphatase were determined using fluorometry as described in previous studies (Bell et al., 2013). The pH, TOC, TN, AP, NO3- - N, NH4+ - N, and SWC were measured using the Qiu et al. (2021) and Xie et al. (2024) techniques. MBC and MBN were determined using the chloroform fumigation extraction method, and their conversion coefficients were 0.45 and 0.54, respectively (Brookes et al., 1985).

Statistical analysis

The AMF communities' α-diversity was assessed using the Quantitative Insights into Microbial Ecology (QIIME, version 2) platform (Zhu et al., 2024). We then averaged the standardized Shannon, Simpson, Chao1, and OTU richness indices to obtain an overall α-diversity index for the AMF communities of A. ordosica and C. korshinskii. The β-diversity of the AMF communities was calculated based on the Bray-Curtis dissimilarity between samples using the ‘adonis’ function in the ‘vegan’ package. The overall differences in AMF communities under each warming treatment were determined using constrained principal coordinate analysis (CPCoA), analysis of similarity (ANOSIM), multiresponse permutation procedures (MRPPs), and permutation multivariate analysis of variance (PERMANOVA), all based on Bray-Curtis dissimilarity. One-way analysis of variance (ANOVA), accompanied by Duncan's post-hoc test, was used to compare the significant differences in soil single function, soil multifunctionality, plant attributes, AMF colonization, and α-diversity under different warming treatments.

Co-occurrence networks for AMF communities under different treatments were constructed using the “WGCNA” package based on the Spearman correlation matrix. Only robust correlations with a Spearman correlation coefficient | ρ | > 0.60 and a p-value < 0.05 after correction by the Benjamini-Hochberg procedure were incorporated into the network analysis. The topological properties of the network were determined using the "igraph" package in R, which included the number of nodes and edges, average clustering coefficient, average path length, network diameter, graph density, network diameter, average degree, and modularity. Finally, the co-occurrence networks of AMF were visualized using Gephi (version 0.9.2) and Cytoscape (version 3.9.1) software.

Using the "Picante" package in R (Kembel et al., 2010), the β-nearest taxon index (βNTI) and Bray–Curtis–based Raup–Crick index (RCbray) were computed to determine the relative importance of deterministic and stochastic processes in AMF community assembly. Briefly, ecological processes were divided into deterministic (homogeneous selection and variable selection) and stochastic (dispersal limitation, homogeneous dispersal, and drift) processes based on βΝΤΙ and RCbray values (Guo et al., 2022). When βNTI < -2 and βNTI > +2 represent community assembly dominated by homogeneous and variable selection, respectively (Stegen et al., 2015). In addition, the relative influences of homogenizing dispersal and dispersal limitation processes were quantified by |βNTI| < 2, RCbray < -0.95, and RCbray > 0.95, respectively, and the influences of “non-dominated” processes were assessed by |βNTI| < 2 and |RCbray| < 0.95 (Stegen et al., 2015).

We performed ordinary least squares linear regressions between AMF community α-diversity, β-diversity, network complexity, assembly processes, mycorrhizal colonization rate, and soil multifunctionality, and conducted Pearson correlation analysis between the α-diversity, β-diversity, network complexity, and single function. In addition, the linear regression correlations between AMF assembly processes, diversity, and network complexity were investigated. We then used partial least squares structural equation modeling (PLS-SEM) to assess the hypothesized direct and indirect relationships between warming degree, AMF diversity, network complexity, community assembly processes, mycorrhizal colonization rate, and soil multifunctionality. The first step in PLS-SEM analysis is to construct an a priori model based on known causal relationships between potential drivers of soil multifunctionality (Hu et al., 2021; Zhu et al., 2024). Before embarking on modeling, it is vital to conduct a bivariate correlation analysis between all the variables to ensure that the relationship between them is sufficiently linear. Finally, the model was evaluated based on fit indices such as Cronbach's alpha, composite reliability, and average variance extracted (AVE), with values greater than 0.7, 0.6, and 0.5, respectively (Zhu et al., 2024), indicating that the model performs well regarding internal consistency, reliability, and convergent validity. In addition, the Goodness of Fit (GoF) index was used to assess the model's overall predictive ability.