Spatial Variation in Soil Fungal Communities across Paddy Fields in Subtropical China

In this work, Illumina MiSeq amplicon sequencing of the ITS region was used to investigate the spatial variation and assembly mechanisms of fungal communities from different soil layers across paddy fields in subtropical China, and the results demonstrate the decreasing importance of environmental filtering and an increase in the importance of dispersal limitation in structuring fungal communities from shallower to deeper soils. Therefore, the results of this study highlight that perceived drivers of fungal community assembly are dependent on sampling depth and suggest that caution is required when interpreting diversity patterns from samples that integrate across depths. This is the first study focusing on assemblages of fungal communities in different soil layers on a relatively large scale, and we thus believe that this study is of great importance to researchers and readers in microbial ecology, especially in microbial biogeography, because the results can provide sampling guidance in future studies of microbial biogeography.

depth-dependent structuring of microbial communities and assemblages. However, rarely has research been conducted to investigate the microbial biogeographical and community patterns of various soil layers across large spatial scales, and certainly not in paddy soils.
If the prevalence of environmental filters, competitive interactions for resources, or barriers to dispersal vary across soil layers, DDRs calculated independently from each layer are likely to be distinct, reflecting changes in the relative importance of environmental factors or more neutral processes shaping the fungal community at different depths. Yet, this is rarely tested. Thus, we investigated the spatial variation and assembly mechanisms of fungal communities from different soil layers across paddy fields in subtropical China. We quantified a range of soil physiochemical properties to reflect environmental heterogeneity and used multiple statistical approaches to disentangle how fungal community assembles in different soil layers. We made the following predictions. (i) The form of DDRs describing how fungal community composition changes with increasing geographic distance between samples is distinct for each of the different soil layers. (ii) The relative importance of environmental factors and spatial factors (reflecting niche-based and more neutral-based mechanisms, respectively) in regulating fungal community assembly changes from shallow to deeper soil layers.

RESULTS
Soil properties. All soil properties except total potassium (TK) exhibited considerable vertical variation (Table 1). Soil organic carbon (SOC), total nitrogen (TN), available nitrogen (AN), carbon-to-nitrogen (C/N) ratio, total phosphorus (TP), available phosphorus (AP), and cation exchange capacity (CEC) decreased significantly with soil depth, while pH and Fe showed the opposite tendency. Soil properties also varied widely across sampling sites in each soil layer (Table 1), while the average coefficient of variation showed no significant difference among soil layers (P Ͼ 0.05; see Table S1A in the supplemental material). Environmental variation (variance-covariance matrix) of standardized soil properties in the 0-to-10-cm (0-10cm)-and 10-20cm-deep layers was significantly higher than that in subsoil (Fig. 1A). Variance-covariance matrix based on all original soil properties confirmed this result (Fig. 1B). To examine how soil properties varied with geographic distance between samples of the same soil depths, we used Mantel correlation tests, and these tests showed that only a few soil properties in the 0-10cm-deep (AP and Fe) and 10-20cm-deep (C/N ratio and Fe) soils covaried with geographical distance. However, many soil properties, including SOC, TN, C/N ratio, AN, Fe, and pH covaried with geographical distance in the 20-40cm-deep soil layer (Table S1A). Additionally, soil parent material did not significantly influence soil properties (Table S1B).
Overall structure of fungal communities. Generally, Ascomycota (46.45%), Zygomycota (31.58%), and Basidiomycota (14.35%) were highly dominant across all samples (see Fig. S2 in the supplemental material). The relative abundance of Ascomycota, Basidiomycota, and Chytridiomycota tended to decrease with soil depth, while Zygomycota showed the opposite trend (Fig. S3). The relative abundance of Glomeromycota and Neocallimastigomycota showed no significant differences across different soil layers (Fig. S3). The relative abundance of Ascomycota and Chytridiomycota was significantly positively correlated with SOC, TN, C/N ratio, and AN but negatively correlated with Fe (P Ͻ 0.05; Table S2). The relative abundance of Zygomycota was negatively correlated with SOC, TN, CN, and AN but positively correlated with Fe (P Ͻ 0.05; Table S2). The relative abundance of Basidiomycota was significantly correlated only with AP (r ϭ 0.387, P Ͻ 0.001) and Fe (r ϭ Ϫ0.415, P Ͻ 0.001), and the relative abundance of Glomeromycota and Neocallimastigomycota was not significantly correlated with any soil properties (Table S2). Fungal biomass and alpha-diversity (richness and Shannon-Wiener index) decreased significantly with soil depth (Fig. 2), whereas spatial variation in fungal biomass and alpha-diversity increased with soil depth (biomass, 27.98%, 27.64%, and 68.15%; average coefficients of richness, 8.73%, 11.03%, and 17.66%; average coefficients of Shannon-Wiener index, 4.81%, 5.67%, and 15.22%; all quoted for 0-10cm-, 10-20cm-, and 20-40cm-deep soil layers, respectively). Spatial variation in fungal community structure (i.e., Sørensen's index, pairwise Bray-Curtis dissimilarity, and Jaccard distance) also increased significantly along with soil depth (Fig. 3A). Principal-coordinate analysis (PCoA) based on beta-diversity indices showed that samples cluster within each soil layer (Fig. S4), and one-way permutational multivariate analysis of variance (per-mANOVA) confirmed that the fungal communities in different soil layers were signifi- cantly different from each other (Table S3A). In addition to soil depth, soil parent material also significantly affected fungal community composition (Table S3B).
Correlation between fungal community and soil variables and geographical distance. Partial Mantel tests showed that fungal beta-diversity indices were significantly correlated with more soil properties in shallower soil layers (0-10cm and 10-20cm), while only a few soil properties showed significant correlations with betadiversity indices in the 20-40cm-deep layer (Table S4). We observed significant correlations between fungal beta-diversity indices and geographical distance between samples in both the 0-10cm and 10-20cm-deep soil layers (Fig. 3B). However, in the 20-40cm layer, no significant DDR could be observed regardless of the beta-diversity index used (Fig. 3B). In addition, between the 0-10cm and 10-20cm soil layers, the DDR slopes showed no significant difference across beta-diversity indices (P Ͼ 0.05).
Contributions of environmental and spatial factors to variation in fungal communities. On the basis of three beta-diversity indices, forward selection procedures were respectively applied to select subsets of environmental and spatial variables which had significant effects on species composition ( Table 2). Indices in the 0-10cm soil layer, regardless of which beta-diversity index was used, the environmental factors pH and TP, and spatial factor PCNM2 (principal coordinate 2 of neighbor matrices) were always selected. In the 10-20cm layer, the environmental factor Fe and the spatial factor PCNM2 could also be frequently selected. In the 20-40cm layer, when using Sørensen's index and Jaccard distance, the environmental factors pH and TN and the spatial factors PCNM2, PCNM11, and PCNM3 were selected. If using Bray-Curtis dissimilarity, no environmental factor could be selected, while the spatial factor PCNM20 had significant effects on these two indices.
Variation partitioning was applied based on the results of forward selection, and in general, the variation in species composition was largely (at least 81.94%) unexplained (Table 3). Of the variance that could be explained, spatial factors rather than environmental factors explained more across all soil layers, except in the 0-10cm soil layer where environmental variables explained slightly more. In the 20-40cm layer, environmental factors barely explain any variance (Table 3).
Community-level habitat niche breadths (Bcom) were estimated, and fungi in the 20-40cm soil layer showed the highest values, followed by those from the 10-20cm layer and finally the 0-10cm layer (Fig. 4C). In contrast to Bcom, the niche overlaps OTUs that occur more frequently than predicted by the model are shown in red, while those that occur less frequently than predicted are shown in yellow. OTUs that occur within prediction are shown in green. Dashed lines represent 95% confidence intervals around the model prediction (black line). (B) The normalized stochasticity test (NST) was developed based on Jaccard distance (NST jac ) with 50% as the boundary point between more deterministic (Ͻ50%) and more stochastic (Ͼ50%) assembly. (C and D) One-way ANOVA and nonparametric Mann-Whitney U test were conducted to test the significance of difference in habitat niche breadth (C) and niche overlaps (D), respectively. among fungi were lowest in the deepest soil layer (20-40cm), while fungi in the shallower layers showed no significant difference in niche overlaps (Fig. 4D).

DISCUSSION
This study quantified spatial variation and drivers of fungal community assembly in paddy field soils from a typical region of subtropical China. Our results consistently revealed that spatial variation in fungi was higher in topsoil than in subsoil and demonstrated the decrease in importance of environmental filtering and an increase in the importance of dispersal limitation in structuring fungal communities from shallower to deeper soils.
We observed obvious spatial variation in fungal communities among our sampling sites (Fig. 3); while the spatial variation was lower than reported from studies conducted on very large scales (34,35), it was higher than that recorded at smaller scales (36,37), suggesting that the degree of spatial variation observed in fungal communities was dependent on the scale. Distance-decay relationships revealed significant correlations between changes in fungal community composition and the geographical distance between samples for 0-10cm and 10-20cm soil layers, but not in the deeper 20-40cm layer (Fig. 3B), suggesting that DDR occurred only in topsoil. Additionally, the slopes of the DDRs from the 0-10cm and 10-20cm soil layers were similar, which may result from the uniform tillage operations that homogenize topsoil to some extent. It should be noted that some factors, such as low environmental heterogeneity, high dispersal, and ecological drift, would greatly weaken the DDR via homogenizing communities (38). Environmental heterogeneity was significantly higher in the shallower soil layers (0-10cm and 10-20cm) than in deeper soil layers, suggesting that low environmental heterogeneity indeed weakened the DDR in the 20-40cm soil layer. Although shallower soil layers (0-10cm and 10-20cm) had higher dispersal rate than deeper soil layers (Fig. 4A), the relatively higher environmental heterogeneity maintained a significant DDR. Judging from DDR and migration rate, we postulate that the fungal communities in shallower soils (0-10cm and 10-20cm deep) are relatively more influenced by environmental factors and dispersal than deeper soils. Additionally, partial Mantel tests showed that a greater number of soil properties were correlated with fungal community composition in shallower soil layers than deeper soil layers (see Table S4 in the supplemental material), whereas few of these soil properties were correlated with geographical distance in these layers (Table S1). This implies that fungal communities in shallower layers are influenced by nonspatially autocorrelated environmental factors.
While higher environmental heterogeneity leads to higher structural heterogeneity of communities (1,14), our study results showed that the degree of environmental variability did not match the extent of community variability ( Fig. 1 and 3A). Thus, environmental factors and environmental variability have less of an impact on variation in species composition compared with other nonenvironmental factors here. Both partial Mantel tests and forward selection demonstrated that some soil properties, especially pH, had significant effects on species composition (Table S4 and Table 2), but their effects were limited and unlikely to be ecologically meaningful. For example, soil pH is always demonstrated to be a key factor affecting microbial assemblages (39), while the coefficients of variation in pH are never higher than 6% in our research, which can keep the soils under relatively strong acidity.
Variation partitioning showed that environmental factors played a slightly more dominant role in driving fungal community assembly in the 0-10cm soil layer. However, spatial factors, rather than environmental factors, played a far larger and more important role in governing fungal community assembly in the deeper soil layers (10-20cm and 20-40cm) ( Table 3), reflecting an increased importance of neutral processes with increasing soil depth. Some other studies on paddy soils found that spatial factors better predicted for fungal community composition compared with environmental factors (5), in contrast to studies of forest soils where environmental factors are generally shown as the better predictors (10). However, our results show that the relative importance of environmental versus spatial factors is soil depth dependent, as environmental factors can better predict fungal communities in the 0-10cm-deep paddy soil layer, while spatial factors are better predictors in deeper soil layers (10-20cm and 20-40cm). Such a pattern should be closely related to the decreased environmental heterogeneity and dispersal rate with soil depth, exposing fungi in topsoil to higher environmental selection and immigration (40). However, a large proportion of variation remains unexplained whichever beta-diversity indices were used (at least 81.94%; Table 3). The unexplained variation may be attributed to unmeasured environmental variables, which could include total dissolved oxygen, the distance from the river or nearby streams, and drainage potential. For example, dissolved oxygen may vary greatly across different soil layers (31,32), promoting redox gradients that may influence the variation in fungal community composition. The distance from sampling site to the nearest river or stream may be another factor affecting fungal community composition, as these bodies of water potentially help the free movement of fungi. In fact, we found it a common phenomenon that a large proportion of variation is rarely explained when using variation partitioning. For example, nearly 90% of variation cannot be explained by spatial or environmental factors on plant community (41), aquatic organisms (42), and soil eukaryotes (33). More importantly, although sampling effects or unmeasured variability may contribute to the unexplained variation, it is tempting to speculate that the high fraction of unexplained variance could be caused by the evolutionary noise produced by ecologically neutral processes of diversification, i.e., through random ecological drift, which cannot be determined by mathematical models (40,43,44).
Dominance test, NST index, and habitat niche breadth were used to help explain the spatial variation in fungal communities and their associated ecological drivers. The dominance test showed that the models had high R 2 values, and more than 80% of species had frequencies within predicted ranges (Fig. 4A), suggesting that the frequency with which fungi occurred in different soil layers can be well described by the neutral model. Even so, some nonneutral process should also be considered. Within each soil layer, there were some fungal species, less than 20%, whose distributions deviated from neutral predictions (Fig. 4A). For example, in the 0-10cm soil layer, the Ͻ20% OTUs occurring outside predictions accounted on average for Ͼ77% of total sequences, while in 10-20cm and 20-40cm soil layers, the Ͻ20% OTUs occurring outside predictions accounted for Ͻ 23% of total sequences on average (Fig. 4A). Random forest models showed that OTUs outside model predictions were more influenced by environmental factors than those inside model predictions (Table S5).These results suggested that more fungi (higher relative abundance, rather than a larger number of taxa) in the 0-10cm layer were selectively enriched or excluded as a result of environmental selection (45). The model also showed a very high migration rate (m) in topsoil and a very low migration rate in subsoil, implying high and unhindered dispersal in topsoil and significant dispersal limitation in subsoil. While some studies suggest that fungi are free to disperse, and thus dispersal limitation does not exist (46), other studies suggest that as fungi are relatively large compared with other microbes (e.g., bacteria) (47), their dispersal may be limited. In our study, the migration rate of fungi was much higher in topsoil. This is easy to monitor, because paddy fields are covered with water for most of the year, and the flow of water would greatly help the free movement of fungi. However, the fungi in subsoil can hardly disperse prior to moving from subsoil to topsoil or, at least, can hardly widely disperse. Additionally, tillage operations that may contribute to the dispersal of fungi are typically only carried out in the topsoil.
The NST index suggested that the relative importance of deterministic processes over stochastic processes in structuring fungal communities decreased with soil depth (Fig. 4B and Fig. S5). The higher environmental heterogeneity in the topsoil exposed soil fungi to a greater range of environmental filters, which drives the unambiguously deterministic process of environmental selection (48). Thus, our results imply that as soil depth increases, environmental selection has an ever lower influence on structuring fungal communities. In contrast, the relative influence of stochastic process in structuring fungal communities increased with soil depth (NST), and this is likely related to increased dispersal limitation (22). Although some stochastic processes such as diversification and ecological drift, which are problematic to quantify, may also increase along with soil depth.
Species with wider niche breadth are considered to be generalists which are less influenced by environmental factors because of higher environmental tolerances (49,50). In our study, habitat niche breadth of fungi continuously increased along with soil depth (Fig. 4C), suggesting that fungi in subsoil with wider niche breadth were governed less by environmental filtering. We initially expected that niche overlaps should be higher in subsoil because of relatively lower resources. However, niche overlaps among fungi were significantly lower in the 20-40cm soil layer than in the 0-10cm and 10-20cm layers (Fig. 4D). This is likely a result of lower fungal biomass and richness (Fig. 2). Fungi in the 20-40cm layer occupied wider niche breadths with lower niche overlaps, suggesting that they can effectively utilize an array of resources with less competition (51) and that they are better adapted to the local environment. Thus, the fungi in the 20-40cm soil layer should be less influenced by deterministic processes, including environmental filtering and competitive exclusion.
Conclusions. We observed obvious spatial variation in fungal communities of paddy fields in subtropical China and found that environmental heterogeneity decreased along with soil depth, while spatial variation in fungal communities showed the opposite tendency. An array of statistical analyses revealed that the fungal community assembly in the 0-10cm-deep layer was primarily governed by environmental filtering and high dispersal, while in the deeper layer (20-40cm), it was primarily governed by dispersal limitation and minimal environmental filtering. Both environmental filtering and dispersal limitation controlled the fungal community assembly in the 10-20cmdeep layer, with dispersal limitation playing the major role. This work highlights that perceived drivers of fungal community assembly are dependent on sampling depth. Thus, future studies interpreting diversity patterns from soil samples that integrate over a wide range of depths should do so with caution, as different ecological mechanisms are likely acting in different soil layers.

MATERIALS AND METHODS
Soil sampling and physicochemical characterization. Soil samples were collected near the end of December 2017 from red paddy soils in Yujiang (Jiangxi Province, China; 116°41= E to 117°09= E, and 28°04= N to 28°37= N), where Ͼ 85% of cultivated land is paddy fields. Sampling sites have subtropical monsoon climates, with abundant sunshine and rainfall (mean annual sunshine hours, 1,739.4 h; mean annual temperature, 17.6°C;mean annual precipitation, 1,750 mm). The total sampling area is 927 km 2 , including 78.2% of hills and 21.2% of plains. The cropping system here is mainly double cropping rice (Oryza sativa L.) (i.e., early and late season rice). Rotary tillage to a depth of 15 to 20 cm (15-20cm) is conducted prior to seedling. The natural conditions, including climate, soil properties, topography, geomorphology, cropping system, and social and economic conditions, including productivity level, are typically representative of subtropical areas of southern China (52).
Sampling sites were chosen to satisfy the following conditions. (i) The whole region needed to be covered. (ii) The main parent material of the soils needed to be included. (iii) Field management, including cropping system and fertilizer applications, should be uniform. On the basis of these principles, 26 sites were selected, with pairwise geographical distances ranging from 1.3 km to 50.7 km (Fig. 5; see Fig. S1 in the supplemental material). The soil samples were collected in December 2017 after the harvest and in the absence of water flooding. Within each site, five 40-cm-deep soil cores (6-cm diameter, free from rice roots) were collected at random locations and partitioned into three depth intervals: 0-10cm, 10-20cm, and 20-40cm. Samples were refrigerated at 4°C using a portable fridge and transported to the laboratory. Samples from each plot were composited by depth, homogenized, and subsampled for subsequent analyses. Subsamples for physical and chemical properties were air dried, ground, and sieved through 2-mm mesh. Subsamples for microbial properties were stored at Ϫ40°C.
Soil chemical properties were determined using the methods described by Pansu and Gautheyrou (53). Soil pH was assayed using a pH meter (FE30; Mettler-Toledo) with 1:2.5 soil-water suspension. Cation exchange capacity (CEC) was determined by saturating the exchange sites of 1 g of each sample twice with 1 M ammonium acetate solution at pH 7, followed by replacing the adsorbed ammonium ions twice with 1 M KCl. Soil organic carbon (SOC) was titrated against 0.5 M ferrous iron solution after it had been digested with 0.8 M K 2 Cr 2 O 4 and concentrated H 2 SO 4 (vol/vol, 1:1) at 150°C for 30 min. Total nitrogen (TN) and available nitrogen (AN) were measured as Kjeldahl N. Briefly, the soil sample was heated and boiled with concentrated H 2 SO 4 . The total nitrogen (TN) was then absorbed by 2% boric acid solution and titrated against 0.1 M sulfuric acid. The available nitrogen was hydrolyzed by 1 M sodium hydroxide and measured by microdiffusion methods. Total phosphorus (TP) and available phosphorus (AP) were extracted with HF-HClO 4 and sodium bicarbonate, respectively, and then determined by the molybdenum blue method using an UV spectrophotometer at 700 nm. Total potassium (TK) was determined using flame emission spectrometry after the soil had been digested in concentrated HF-HClO 4 (vol/vol, 2:1). Free iron (Fe) of the soil was extracted by dithionite-citrate-bicarbonate (DCB) solution with 1 g of soil being digested in 40 ml of 0.3 M sodium citrate and 5 ml of 1 M sodium hydrogen carbonate at 353 K for 30 min, and the amount of free iron was then determined by flame atomic absorption spectrophotometry.
Soil DNA extraction, amplification, Illumina sequencing, and sequence processing. Soil DNA was extracted from 0.5 g of soil (fresh weight) using a FastDNA SPIN kit (MP Biomedicals, CA, USA) and then subsequently purified using a PowerClean DNA clean-up kit (MoBio, CA, USA) according to the manufacturers' instructions. The concentration and quality of the extracted DNA were measured using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, DE, USA). Quantitative PCR was done on the Bio-Rad CFX96 touch real-time PCR detection system following protocols previously described (54). Results are reported as gene copy numbers and are used to reflect fungal biomass.
Each of the 78 DNA samples was amplified separately using the fungal PCR primers ITS1F (5=-CTTG GTCATTTAGAGGAAGTAA-3=) and ITS2 (5=-GCTGCGTTCTTCATCGATGC-3=) (55) that target the internal transcribed spacer 1 (ITS1) region. PCR products were then sequenced on the Illumina MiSeq PE250 platform. Raw sequence data were analyzed using the Quantitative Insights into Microbial Ecology (QIIME) pipeline (v1.9.1) (http://qiime.org/) (56). Paired-end reads were merged using FLASH (57). Reads with length of Ͻ200 bp or with average quality scores of Ͻ25 were removed, resulting in 796,863 high-quality sequences. ITSx 1.0.11 (http://microbiology.se/software/) was then used to remove 5.8S and 28S regions from merged sequences (58). Any chimeric sequences were removed using the USEARCH tool based on the UCHIME algorithm (59). Operational taxonomic unit (OTU) picking was performed using the "pick_otus.py" command with the nondefault UCLUST algorithm (the parameters were as follows: picking method, uclust; similarity of 0.97; max_accepts of 20, and max_rejects of 100) (60). Sequences were clustered into 12,474 OTUs after excluding singletons and rarefying to 44,296 sequences per sample (based on the sample with the minimum numbers of reads) (61,62). The taxonomic identity of each OTU was then determined based on comparisons against the UNITE database (v7) (https://unite .ut.ee/). Statistical analysis. Alpha-diversity indices, including richness and Shannon-Wiener index were calculated in QIIME using the "alpha_diversity.py" script. Statistically significant differences in soil properties, fungal biomass, and alpha-diversity indices were determined by one-way analysis of variance (ANOVA), along with the use of Duncan's test for multiple comparisons (P Ͻ 0.05). If the variances of observations were heterogeneous, nonparametric Mann-Whitney U test were used to determine the statistical significance. Variance-covariance matrix based on all soil properties was calculated to indicate the overall variation in soil properties. Three beta-diversity indices, including Sørensen's index, pairwise Bray-Curtis dissimilarity, and Jaccard distance coupled with principal-coordinate analysis (PCoA) were conducted to indicate the community dissimilarities. Partial Mantel tests were conducted to determine the potential effects of each soil property on fungal composition. Variation-partitioning analysis was conducted to disentangle the relative importance of environmental factors and spatial factors on variation in fungal communities. Spatial variables were derived from the principal coordinates of neighbor matrices (PCNM) algorithm, which was able to deconvolute total spatial variation into a discrete set of explanatory spatial scales (63). Forward selection procedures were subsequently used to select respective subsets of environmental and spatial variables. The forward selection was stopped if the significance level (P Ͼ 0.05) was reached, or if no improvement of selection criterion (R 2 ) was seen when adding any additional variables. A two-way permutational multivariate analysis of variance (permANOVA) was then performed with the selected variables using the R script provided by Wu et al. (38). Pure environmental variation without a spatial component represents the strength of environmental filtering, while pure spatial variation without an environmental component is interpreted as the effect of dispersal limitation. The fractions of explained variance are based on adjusted fractions (R 2 adj , adjusted coefficient of multiple determination), which accounts for the number of variables and sample sizes. The significance of each component via partitioning was evaluated with a permutation test, except for the interaction term and residuals (these cannot be tested statistically).
A neutral assembly model (the so-called dominance test) was used to determine the potential contribution of neutral processes to the community assembly by predicting the relationship between the occurrence frequency of OTUs and their relative abundance (64). This model evaluates whether the microbial assembly process from a metacommunity follows a neutral model (inside model predictions) or a niche-based process (outside model predictions) as a function of the metacommunity log abundance. Random forest analyses were subsequently performed to quantitatively evaluate the importance of predictors influencing OTUs that occurred outside or inside predictions of the dominance test. The importance of each predictor was determined by assessing the decrease in prediction accuracy (that is, the increase in the mean square error [MSE] between observations and predictions) when the data for the predictor were randomly permuted. This decrease was averaged over all trees to produce the final measure of importance. These analyses were conducted using the "randomForest" package of the R statistical language (65). The significance of predictor importance was assessed by using the "rfPermute" package.
We further applied the normalized stochasticity ratio (NST) to help confirm fungal community assembly processes. NST is an index developed with 50% as the boundary point between more deterministic (Ͻ50%) and more stochastic (Ͼ50%) assembly (66). We choose NST to indicate assembly processes because our research met the requirements of this method: (i) local/landscape scale sampling as opposed to global scale; (ii) n Ն 6. This analysis was conducted in the R statistical language (65) using "NST" package (the parameters were as follows: "dist.method" of "bray"/"jaccard," "abundance.weighted" of "TRUE", and "rand" of "1000"). By considering the overall performance of similarity metrics, NST based on Jaccard distance (NST jac ) is recommended for estimating the magnitude of stochasticity in community assembly (66), but NST based on Bray-Curtis dissimilarity (NST bray ) is also calculated in our research to further verify NST jac .
Niche breadth and niche overlaps were respectively calculated according to Levin's niche breadth index and Levin's niche overlap index (38). Briefly, Levin's niche breadth index was determined as follows: where B j represents the habitat niche breadth of OTU j in a metacommunity, N is the total number of communities of each metacommunity, and P ij is the proportion of OTU j in community i. A high B indicates that the OTU occurs widely and evenly along a wide range of locations, representing wide habitat niche breadth. We calculated the average B values from all taxa in a single community (Bcom) as an indicator of habitat niche breadth at the community level. Levins' niche overlap index (O) was calculated as follows: where O jk represents the niche overlap between OTU j and OTU k , N is the total number of communities of each metacommunity, P ij is the proportion of OTU j in community i, P kj is the proportion of OTU k in community i. A high O indicates that the species exhibited more niche overlap.
Data accessibility. The ITS sequences used in this study were submitted to the NCBI Sequence Read Archive (SRA) under the accession number SRP200912.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.