Decreased soil multifunctionality is associated with altered microbial network properties under precipitation reduction in a semiarid grassland

Our results reveal different responses of soil multifunctionality to increased and decreased precipitation. By linking microbial network properties to soil functions, we also show that network complexity and potentially competitive interactions are key drivers of soil multifunctionality.


Appendix 1 Method for the determination of soil organic carbon fraction
POXC was determined as described by Blair et al. [1].Air-dried soil (containing 15 mg of C) was added to tube, and shaken with 20 ml of 333 mM potassium permanganate for 1 h at 180 rpm.The tube was then centrifuged for 10 min at 4,000 g and the supernatant solutions diluted by 1:250 with distilled water.The absorbance of diluted solutions was measured at 565 nm.HCl-ROC was determined as described by [2]. 1 g air-dried soil was digested with 6 M HCl at 125°C for 16 hours.The residue after the hydrolysis was separated from the supernatant by centrifugation (10 min at 6,037 g), and washed five times with de-deionized water to remove any chlorine residue.
The remaining residue was oven-dried at 60 °C, finally the carbon content is measured by the K2Cr2O7 oxidation method.
The labile and recalcitrant carbon functional groups of soil C were determined by the following methods as described by Rui et al. [3].Samples were finely ground, mixed with potassium bromide (IR grade; soil: KBr ratio = 1:100, w:w), and analyzed using a Thermo Nicolet 6700 infrared spectrometer (Thermo Electron Scientific Instruments Corp).Reflectance spectra (4,000 -400 cm −1 ) were collected with 64 scans and a resolution of 4 cm −1 .We evaluated relative peak areas for the following organic functional groups: aliphatic C-H functional group of methyl and methylene groups (ALC: wavenumber 3,010 to 2,800 cm −1 ), C-O in both poly-alcoholic and ether functional groups (PEC: wavenumber 1,170 to 1,148 cm −1 ), aromatic C = C stretch and/or asymmetric -COO-stretch (AC: wavenumber 1,660 to 1,580 cm −1 ), carboxylic C=O (CC: wavenumber 1,740 to 1,700 cm −1 ).Given that carbohydrate and aliphatic functional groups are assimilated more easily than aromatic and carboxylic compounds [4].we categorized ALC and PEC as labile carbon structures, whereas AC and CC were categorized as recalcitrant carbon structures.

A-2.1 Nutrient provisioning
Soil dissolved organic carbon (DOC) and nitrogen (DON) were extracted with MilliQ water and measured using a TOC analyzer (Shimadzu, Kyoto, Japan).Soil inorganic N (NH4 + -N and NO3 --N) and available phosphorus (AP) was extracted with 2 M KCl solution and 0.5 M NaHCO3, respectively, and then were measured using a continuous flow analyser (AutoAnalyzer 3 SEAL, Bran and Luebbe Corp).

A-2.2 OM decomposition
The activities of the five extracellular enzymes (i.e., BG, CBH, NAG, LAP, and ALP) were measured fluorometrically using a 200-μM solution of substrates labelled with 4-methylumbelliferone (4-MUB) or 7-amino-4-methylcoumarin (7-AMC) [5].To this end, we used a 96-well fluorometric microplate, which has 16 analytical replicates for each soil sample and 8 analytical replicates for blanks, negative controls, reference standards, and quench standards.Specifically, one gram of fresh soil sample was homogenized in 125 ml sodium acetate buffer (pH = 8.The activities of soil oxidases (PER and PO) were measured by the colorimetrical method using L-dihydroxyphenylalanine (L-DOPA) and L-DOPA + H2O2 as substrates, respectively [6].Briefly, 1 g of frozen soil was homogeneously dispersed in 125 ml of sodium acetate buffer (pH = 8.5) on a constant-temperature shaker (25 °C) at 180 rpm for 2 h.Eight aliquots (150 μl each) of the prepared soil suspension and L-DOPA solution (50 μl each, 5 mM) were successively added into a clear 96-well microplate using an eight-channel micropipette.For PER, another 10 μl of 0.3% H2O2 were dispensed into the wells combined with the suspensions and substrates.The microplate was incubated in the dark at 25 °C for 5 h, and then 10 μl of 1 M NaOH was added to stop the reaction.Absorbance at 450 nm was measured using a microplate spectrophotometer (Tecan Infinite M200, Salzburg, Austria).

A-2.3 Microbial growth efficiency
Soil microbial biomass C (MBC), N (MBN), and P (MBP) were measured using the chloroform fumigation-extraction techniques [7−9].We calculated the microbial carbon use efficiency (CUE) based on the extracellular enzyme stoichiometric model [10], and further incorporated soil respiration data to calculate biomass turnover rate (BTR).Microorganisms with higher CUE can transformed more C into biomass and stabilized sequester in soils [11].Moreover, BTR also determines the fate of C, with higher turnover rates generally promoting biomass accumulation [12].These two microbial parameters underlie critical mechanisms in the global C cycle, especially in a climate-changing world [13,14].We performed the calculation by the following three equations [10]. (1) Where

A-2.4 Soil respiration
To measure soil respiration (Rs), in one subplot of each plot, a thin-walled polyvinyl chloride (PVC) respiration collar (20 cm diameter and 5 cm height) was permanently inserted approximately 3 cm into the soil to minimize the potential mechanical disturbance during the experiment establishment.Soil respiration were recorded using an IRGA infrared gas analyser (Li-8100; LiCor) monthly during May-September in 2019.While measuring soil respiration, we also measured soil temperature (ST) and water content (SWC) at the 10 cm depth using a portable temperature meter and a TDR-100 (Spectrum Technologies).Each measurement is selected to start at 9.00-11.00am on a clear day and last for three days.We measured 18 plots in the same sequence to reduce temporal and spatial biases.In addition, we frequently clipped living plants inside the PVC collars and left the clipped plant materials on the soil within the collars.Each clipping was finished at least 1 day before the soil respiration measurement to minimize disturbance influence.We used the mean (n=15) of all soil respiration data during the measurement period to represent soil respiration during growing season in 2019.

A-2.5 Soil multifunctionality
As a human construct and multidimensional measure of ecosystem functions, multifunctionality quantifies the ability of an ecosystem to simultaneously deliver multiple ecosystem processes and services [15].In this study, we evaluated multifunctionality through multiple approaches.
(1) multiple single function [16]: multiple single functions represent the actual magnitude of multifunctionality in each metric.It is derived from the value after standardization (a common scale ranging from 0 to 1) of each directly measured metric; (2) the averaging multifunctionality index [17]: this index is obtained by averaging the standardized values of individual functions and is the most widely used and easily available multifunctionality index; (3) weighted multifunctionality by four groups of ecosystem function [18]: this index averages ecosystem functions into four ecosystem services: nutrient provisioning, microbial growth efficiency, LOM decomposition, and ROM decomposition before calculating multifunctionality, so that each ecosystem service contributes equally to multifunctionality; (4) the principal coordinate multifunctionality index [15]: this index characterizes the different dimensions of multifunctionality by extracting the principal axes from the principal coordinate analysis of all the single functions.further use [19].To profile soil bacterial communities, we amplified the V4 hypervariable region of 16S rRNA gene with the primer sets 515F and 907R [20].For fungal communities, we amplified the ITS1 region with the primers sets ITS1-ITS2 [21].Amplification was performed with an ABI GeneAmp® 9700 PCR thermocycler (ABI, CA, USA).
To permit multiplexing of samples, a 10-bp barcode unique to each sample was attached to the 5′ end of primers.The reaction system of bacterial PCR amplification contained 0.4 μl of each primer, 1.25 μl of template DNA (10 ng), 0.4 μl of FastPfu polymerase (Beijing TransGen Biotech Co., Ltd, China), 2.5 μl of 10 × PCR buffer, 2.0 μl of dNTPs (2.5 mM), 5.0 μl of 5 × high enhancer, and 13.45 μl of sterile ultrapure water, in a total reaction volume of 25 μl.The PCR amplification of 16S rRNA was performed at an initial denaturation temperature of 95 °C for 3 min, followed by 27 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 45 s.A final elongation step at 72 °C was run for 10 min.The PCR reaction system of the ITS gene was also performed in a 25 μl mixture, which contained 0.5 μl of the two primers (30 μmol μl −1 ), 1.5 μl of template DNA (10 ng), and 22.5 μl of Platinum PCR SuperMix (Invitrogen, Shanghai, China).The thermal cycling program consisted of an initial annealing temperature at 95 °C for 2 min, followed by 30 cycles with denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, extension at 72 °C for 30 s, and final extension at 72 °C for 5 min.
The three replicated amplification products of bacteria and fungi in each soil sample were mixed to provide one final PCR product.Each mixed PCR product was purified using the Qiagen™ PCR purification kit (Qiagen, Valencia, California, United States) following the manufacturer's protocol, and then were eluted in sterile water.
Subsequently, 2.0% agarose gel electrophoresis were used to quantified the concentration of each mixed gene (16S rRNA gene and ITS rRNA gene).Finally, the purified amplicons were pooled in equimolar and paired-end sequenced (2 × 300) on an Illumina MiSeq platform (Majorbio, Shanghai, China) according to the standard protocols.

A-3.2 Sequence processing
Raw FASTQ files of sequencing were de-multiplexed and quality-filtered using Quantitative Insights Into Microbial Ecology (QIIME, Version 1.9.1) workflow [19].
Briefly, the tags sequence, primer sequence, the sequences < 200 bp, an average quality score of < 25, and the reads containing ambiguous bases or any unresolved nucleotides were removed to obtain the high-quality sequence for subsequent analysis.Forward and reverse reads were incorporated into full-length sequences with FLASH (Fast Length Adjustment of Short reads).The above processed sequences were clustered into operational taxonomic units (OTUs) defined by 97% similarity using the UPARSE 11 pipeline [22], and chimeras were eliminated during this procedure.Subsequently, the largest number of sequences in each OTU were selected as the representative sequence to prepare the OTUs table.The ribosomal database project (RDP 11.

A-4.1 Network construction
All microbial co-occurrence networks were constructed on the basis of Pearson correlations of log-transformed OTU abundances, and then an RMT-based approach was used to automatically determine the cut-off threshold for the correlations.The RMT-based network approach has several advantages [23,24].First, it has a sound theoretical foundation, which is based on two general laws of random matrix theory.
Second, it has a mathematically defined non-arbitrary correlation cutoff, which minimizes the uncertainty in network construction and comparison.Third, it has a high tolerance for noise from non-random, system-specific features.Moreover, it provides broad applicability to ecological datasets in different formats.Therefore, we chose to build the network based on this approach to minimize the interference factors in network construction and comparison.
To avoid potential spurious associations of rare OTUs from affecting reliability, data filtering was performed prior to correlation calculation.For the metacommunity cross-kingdom co-occurrence network, only OTUs present in 9 of the 18 samples were included for calculation.For the independent networks under each precipitation regime, only OTUs present in 4 of the 6 samples were included.The Pearson correlation coefficient was calculated based on the log-transformed relative abundances of OTUs.
The correlation matrix obtained from the analysis of the RMT-based network approach was used to determine the correlation threshold for network construction.Only robust correlations with absolute values greater than or equal to the threshold were used to construct the network.
To test whether the derived empirical network was non-random and scale-free, the networks were evaluated against 100 random networks generated by randomly rewiring the links among the nodes [25].With each randomization, the same suite of network topological properties was calculated.The means and standard deviations of these properties from the 100 randomizations were calculated and compared with those from the corresponding empirical MENs.Network randomization was performed in the Molecular Ecological Network Analysis Pipeline.

A-4.2 Network characterization
For the metacommunity cross-kingdom co-occurrence network, we extracted subnetworks by preserving the phylotypes (OTU) of individual soil samples using the induced_subgraph function in "igraph" package in R.Then, to describe and compare the topology of each subnetwork, we calculate a set of metrics: the number of nodes and edges, average degree, graph density, average path length, network diameter, clustering coefficient, and betweenness centrality.Each node represents one OTU, and each edge represents a strong and significant correlation between two nodes; average degree refers to the average connections of each node with another unique node in the network; graph density refers to the intensity of connections among nodes; Average path length refers to the average network distance between all pairs of nodes; clustering coefficient represents the degree to which the nodes tend to cluster together; and graph density refers to the intensity of connections among nodes betweenness centrality refers to the potential influence of a particular node on the connections of other nodes.
Therefore, higher numbers of nodes and edges, average degree, clustering coefficient, betweenness centrality, and graph density, and lower average path lengths suggest a more connected network, reflecting more potential complexity of soil networks.Finally, these topological features of the soil meta-community network are used to calculate the network complexity of each sample.
For the independent networks under each precipitation regime, we calculated network robustness and vulnerability to understand the effects of altered precipitation on network stability.The robustness (the resistance to node loss) of a given network is defined as proportion of the remaining species in this network after random removal to simulate species extinction.To test the effects of species removal on the remaining species, we calculated the abundance-weighted mean interaction strength (wMIS) of node i as In the equation, where bj is the relative abundance of species j and sij is the association strength between species i and j, which is measured by Pearson correlation coefficient.The proportion of the remaining nodes was reported as the network robustness.The robustness value was defined when 50% of random nodes were removed.
The vulnerability of each node measures the relative contribution of the node to the global efficiency.The vulnerability of a network is indicated by the maximal vulnerability of nodes in the network as Where E is the global efficiency and Ei is the global efficiency after removing node i and its entire links.The global efficiency of a graph was calculated as the average of the efficiencies over all pairs of nodes: Where d (i,j) is the number of edges in the shortest path between node i and j.
We also tested whether the networks were scale-free and modular.According to classification criteria [26], all networks exhibited scale-free properties (R 2 of powerlaw degree distribution ranging from 0.754 to 0.861).Modularity (M) measures the degree to which a network is compartmentalized into different modules, in which the nodes within a module are highly connected but have very few connections with nodes from other modules.Since the size and connectivity of different networks vary greatly, relative modularity (RM: a measure of how modular a network is compared to the average expected modularity) is more meaningful for comparing the modular structures of different networks [27].In this study, RM was calculated as the ratio of the difference between the modularity of an empirical network and the mean of modularity from the random networks over the mean of modularity from the random networks: where M is the modularity of the empirical network, and  is the mean of modularity from the random networks (Table S1).All the networks were highly modular based on both the modularity (> 0.4) and RM (> 0, meaning that modularities of empirical networks were larger than those of random networks).

This supplementary file includes: 1 .
Figure S1 Changes in the annual precipitation on the Loess Plateau from 1901 to 2100.The data used for the fitting came from http://gdo-dcp.ucllnl.org/downscaled_cMIp_projections/.

Figure
Figure S2 Rarefaction curves for the richness of bacteria and fungi.Each line represents a soil sample.

Figure
Figure S3 Comparison of organic carbon fractions under different treatments.Vertical bars denote standard errors of mean values (n = 6).Different letters above the vertical bars indicate significant differences among treatments at P < 0.05 (multiple comparison with Kruskal-Wallis tests).DP, 50% decrease in mean annual precipitation; Control, ambient precipitation; IP, 50% increase in mean annual precipitation.

Figure
Figure S4 Correlations between multiple functional metrics.(A-B) Linear relationships between averaging soil multifunctionality and weighted multifunctionality as well as the multifunctionality represented by the first axes from a principal coordinate analysis including the 17 functions.(C) Spearman correlations between the single soil functions and the different components (Comp) of multifunctionality represented by the multiple axes from a principal coordinate analysis.Significance levels were as follows: *P < 0.05, **P < 0.01.

Figure
Figure S5 Relationships between soil function and microbial diversity using various measures.(A) Relationships estimated by Spearman correlation.(B) Functional effects (standardized slopes) of bacteria (blue circle), fungi (black square) and their average

Figure
Figure S6 Relationship between single soil function and microbial diversity.(A) Spearman correlations between single soil function and microbial biodiversity.(B) Functional effects (standardized slopes) of bacteria (blue circle), fungi (black square) and their average diversity (red triangle) index on single soil functions.Significance levels were as follows: *P < 0.05, **P < 0.01.

Figure
Figure S7 Relative importance of factors for multiple soil function indices as indicated by a random forest model.Percentage increases in the MSE (mean squared error) of variables was used to estimate the importance of these predictors, and higher MSE% values implied more important predictors.Significance levels were as follows: *P < 0.05, **P < 0.01.

Figure
Figure S8 Normalized stochasticity ratio and neutral model fitting of bacterial and fungal community assembly.(A) The magnitude of normalized stochasticity ratio (NST) and mean habitat niche width of bacteria and fungi.The value of 50% as the boundary point between more deterministic (< 50%) and more stochastic (> 50%) assembly.Vertical bars denote standard errors of mean values.Different letters above the vertical bars indicate significant differences among treatments at P < 0.05 (multiple comparison with Kruskal-Wallis tests).(c) fit of the neutral model in microbial communities across different treatment.m indicates the estimated migration rate and R 2 indicates the fit to the neutral model.DP, 50% decrease in mean annual precipitation; Control, ambient precipitation; IP, 50% increase in mean annual precipitation.

Figure
Figure S9 Spearman correlations between single soil functions and the diversity (number of OTUs) of soil phylotypes within the main network clusters.Significance levels were as follows: *P < 0.05, **P < 0.01.

Figure
Figure S10 Total OTU richness for all samples at the phyla level within the major network clusters.

Figure
Figure S11 Principal coordinate analysis of network topologies with different treatments.

Figure
Figure S12 Bacterial co-occurrence network (A) and features (B) under different precipitation regimes.Error bar corresponds to the standard deviation of 100 repetitions of the simulation.Different letters above the vertical bars indicate significant differences among treatments at P < 0.05 (multiple comparison with Kruskal-Wallis tests).DP, 50% decrease in mean annual precipitation; Control, ambient precipitation; IP, 50% increase in mean annual precipitation.

Figure
Figure S13 Fungal co-occurrence network (A) and features (B) under different precipitation regimes.Error bar corresponds to the standard deviation of 100 repetitions of the simulation.Different letters above the vertical bars indicate significant differences among treatments at P < 0.05 (multiple comparison with Kruskal-Wallis tests).DP, 50% decrease in mean annual precipitation; Control, ambient precipitation; IP, 50% increase in mean annual precipitation.

Figure
Figure S14 Cross-kingdom co-occurrence network (A) and features (B) under different precipitation regimes.Error bar corresponds to the standard deviation of 100 repetitions of the simulation.Different letters above the vertical bars indicate significant differences among treatments at P < 0.05 (multiple comparison with Kruskal-Wallis tests).DP, 50% decrease in mean annual precipitation; Control, ambient precipitation; IP, 50% increase in mean annual precipitation.

Figure
Figure S15 Standardizing direct, indirect, and total effects of different factors on soil multifunctionality from structural equation model (SEM).
5) on a constant-temperature shaker (25 °C) at 180 rpm for 2 h.And then, 200 μl soil slurry and 50 μl of 200 μM fluorometric substrate were added to the microplate wells that served as the sample assay.In addition, 200 μl buffer + 50 μl of 200 μM fluorometric substrate, 200 μl buffer + 50 μl of 10 μM fluorescence standard solution, 200 μl soil slurry + 50 μl buffer, and 200 μl soil slurry + 50 μl of 10 μM fluorescence standard solution were also dispensed into the wells that served as the negative controls, reference standards, blanks, and quench standards, respectively.The microplates were incubated in the dark at 25 °C for 4 h.After incubation, 10 μl of 0.5 M NaOH was added to each well to terminate the reactions, and the fluorescence values were measured using a microplate fluorometer (Tecan Infinite M200, Salzburg, Austria) with the excitation and emission filters of 365 and 450 nm, respectively.
CUEC:X is the CUE calculated based on C:N or C:P ratio; DOC/AX is DOC:DON or DOC:AP; KX is the half-saturation constant with a value of 0.5; CUEmax represents the upper limit for microbial growth efficiency (0.6) based on thermodynamic constraints.qCO2 (d −1 ) is the microbial respiration rate per unit biomass.Note that although extracellular enzyme-based stoichiometric models are widely used to estimate microbial utilization efficiency, they have inherent disadvantages compared to isotope-based methods.Therefore, subsequent studies are encouraged to use isotope-based methods to accurately quantify microbial carbon use efficiency.

Appendix 3
Method for DNA extraction and MiSeq sequencing A-3.1 Soil DNA extraction, PCR amplification, and Illumina sequencing Soil DNA was extracted from 0.5 g fresh soil of each sample using a FastDNA SPIN Kit for Soil (MP Biochemicals, Solon, OH, USA).The quality and concentration of DNA were determined by 1.0% agarose gel electrophoresis and a NanoDrop® ND-2000 spectrophotometer (Thermo Scientific Inc., USA) and kept at -80 ℃ prior to 5) classifier was used to assign obtained bacterial and fungal OTUs to taxonomic groups based on the SILVA database (version 138) and the UNITE fungal ITS database (version 8.0), respectively.The complete dataset was sent to the Sequence Read Archive (SRA) database of the NCBI under the accession numbers of PRJNA870174 for bacteria and PRJNA870219 for fungi.A-3.3Data processingAfter quality sequencing, both bacterial diversity (959,199 high-quality sequences) and fungal diversity (1,068,716 high-quality sequences) were obtained with the 515F/907R (bacterial 16S rRNA) and ITS1-ITS2 (fungal ITS) primer sets across all soil samples, respectively.The number of bacterial sequences varied from 39,702 to 58,298 per sample, whereas the number of fungal sequences varied from 46,495 to 66,192 per sample.To minimize any bias in the distribution of taxa, the bacterial and fungal diversities of each treatment were calculated based on randomly selected sequences until the count reached saturation in the rarefaction curves.For the downstream analysis of soil bacteria, datasets were rarefied to 39,702 sequences, whereas for the downstream analysis of soil fungi, datasets were rarefied to 46,495 sequences.The OTUs were clustered at a similarity above 97% to calculate the diversity indices.
Figure S1 Changes in the annual precipitation on the Loess Plateau from 1901 to 2100.The data used for the fitting came from http://gdo-dcp.ucllnl.org/downscaled_cMIp_projections/.

Figure S2
Figure S2 Rarefaction curves for the richness of bacteria and fungi.Each line represents a soil sample.

Figure S3
Figure S3 Comparison of organic carbon fractions under different treatments.Vertical bars denote standard errors of mean values (n = 6).Different letters above the vertical bars indicate significant differences among treatments at P < 0.05 (multiple comparison with Kruskal-Wallis tests).DP, 50% decrease in mean annual precipitation; Control, ambient precipitation; IP, 50% increase in mean annual precipitation.

Figure S4
Figure S4 Correlations between multiple functional metrics.(A-B) Linear relationships between averaging soil multifunctionality and weighted multifunctionality as well as the multifunctionality represented by the first axes from a principal coordinate analysis including the 17 functions.(C) Spearman correlations between the single soil functions and the different components (Comp) of multifunctionality represented by the multiple axes from a principal coordinate analysis.Significance levels were as follows: *P < 0.05, **P < 0.01.

Figure
Figure S5 Relationships between soil function and microbial diversity using various measures.(A) Relationships estimated by Spearman correlation.(B) Functional effects (standardized slopes) of bacteria (blue circle), fungi (black square) and their average diversity (red triangle) index on soil function.Significance levels were as follows: *P < 0.05, **P < 0.01.

Figure S6
Figure S6 Relationship between single soil function and microbial diversity.(A) Spearman correlations between single soil function and microbial biodiversity.(B) Functional effects (standardized slopes) of bacteria (blue circle), fungi (black square) and their average diversity (red triangle) index on single soil functions.Significance levels were as follows: *P < 0.05, **P < 0.01.

Figure
Figure S7 Relative importance of factors for multiple soil function indices as indicated by a random forest model.Percentage increases in the MSE (mean squared error) of variables was used to estimate the importance of these predictors, and higher MSE% values implied more important predictors.Significance levels were as follows: *P < 0.05, **P < 0.01.

Figure S8
Figure S8 Normalized stochasticity ratio (A), mean habitat niche width (B), and neutral model fitting (C) of bacterial and fungal community assembly.The value of 50% as the boundary point between more deterministic (< 50%) and more stochastic (> 50%) assembly.Vertical bars denote standard errors of mean values.Different letters above the vertical bars indicate significant differences among treatments at P < 0.05 (multiple comparison with Kruskal-Wallis tests).In Panel C, m indicates the estimated migration rate and R 2 indicates the fit to the neutral model.DP, 50% decrease in mean annual precipitation; Control, ambient precipitation; IP, 50% increase in mean annual precipitation.

Figure S9
Figure S9 Spearman correlations between single soil functions and the diversity (number of OTUs) of soil phylotypes within the main network clusters.Significance levels were as follows: *P < 0.05, **P < 0.01.

Figure
Figure S10 Total OTU richness for all samples at the phyla level within the major network clusters.

Figure S11 Principal coordinate analysis of network topologies with different treatments . Figure S12
Figure S11 Principal coordinate analysis of network topologies with different treatments.

Figure
Figure S13 Fungal co-occurrence network (A) and features (B) under different precipitation regimes.Error bar corresponds to the standard deviation of 100 repetitions of the simulation.Different letters above the vertical bars indicate significant differences among treatments at P < 0.05 (multiple comparison with Kruskal-Wallis tests).DP, 50% decrease in mean annual precipitation; Control, ambient precipitation; IP, 50% increase in mean annual precipitation.

Figure
Figure S14 Cross-kingdom co-occurrence network (A) and features (B) under different precipitation regimes.Error bar corresponds to the standard deviation of 100 repetitions of the simulation.Different letters above the vertical bars indicate significant differences among treatments at P < 0.05 (multiple comparison with Kruskal-Wallis tests).DP, 50% decrease in mean annual precipitation; Control, ambient precipitation; IP, 50% increase in mean annual precipitation.

Figure
Figure S15 Standardizing direct, indirect, and total effects of different factors on soil multifunctionality from structural equation model (SEM).