The experimental area is located in northwest China between 95°1′ and 108°18′E latitude and 36°40′ and 40°19′N longitude and is the predominant arid area in China. The mean annual temperature in the region ranges from 1.6 to 7°C, and the mean annual precipitation ranges from 50 to 250 mm (Table 1).

Three arid regions were selected from different provinces, including Haixi in Qinghai Province (HX), Jiuquan in Gansu Province (JQ) and Etoke in Inner Mongolia Province (ET) in northwest China (Figure 1A). The climates are classified as continental plateau climate, continental arid climate and temperate arid climate for HX, JQ and ET, respectively. The ecosystem types are alpine desert grassland, desert grassland and desert grassland for HX, JQ and ET, respectively. The soil types are alpine meadow soil, aeolian sand soil and meadow chestnut soil for HX, JQ and ET, respectively. In each study region (Figure 1B), 24 soil plots were selected (only 21 plots in JQ due to the loss of 3 samples) and all plots were divided into two habitats based on the dominant plant species: Artemisia habitat (37 plots) dominated by Artemisia ordosica and Artemisia annua and Poaceae habitat (32 plots) dominated by Stipa spp., Leymus spp., and Achnatherum spp. The distance between sites in a region is approximately 20 km. At each plot 5 soil cores (5 cm diametre) from the upper 10 cm were sampled and combined into a single bulk sample, and a total of 69 bulk samples were collected for further analysis. All soil samples were sampled in August 2020 and divided into two subsamples and transported to the laboratory on ice. One subsample was stored at 4°C for measuring soil physicochemical properties, and the other was kept at −80°C for measuring the soil microbial matrix.

The air-dried soils were used to measure soil total organic matter (SOM) with K₂Cr₂O7–H₂SO4 titrimetric method, soil total nitrogen (TN) and total phosphorus (TP) contents using Kjeldahl digestion and vanadium molybdate yellow colorimetry, respectively. Fresh soils (2.5 g) were extracted with 2 M KCl and filtered to determine soil nitrate (NO₃⁻) and ammonium (NH₄⁺) contents using Clever Chem 200+. The particle size composition was determined using a laser particle size analyser (mastersizer2000, Malvern Instruments, United Kingdom). Soil electrical conductivity was measured with a DDS-307 conductivity metre (INESA Scientific Instrument, China). Soil pH was determined using a Sartorius pH meter (PB-10, from Sartorius Corporate Administration GmbH in Göttingen, Germany) by mixing fresh soil with water in a 1:5 ratio. The soil soluble salt content (Salinity) was measured using the residue drying method (Zhang et al., 2019). Soil water content (SWC) was determined as the ratio of the mass of water in the soil sample to the mass of dried soil. Field water capacity (FWC) was calculated as the maximum water content divided by the volume of the soil core. Soil bulk density (SBD) was measured by dividing the mass of oven-dried soil by the volume of the soil core. The saturated hydraulic conductivity (KS) was estimated in the laboratory using a permeameter (Eijkelkamp Agrisearch Equipments, The Netherlands), and the constant head method was used. Data for mean annual air temperature (MAT) and mean annual precipitation (MAP) were obtained from the WorldClim database, while the potential evapotranspiration was obtained from a monthly 1 km dataset (Shouzhang, 2022). The arid index was calculated by dividing the potential evapotranspiration by the mean annual air temperature. The plant index was calculated using Normalized Difference Vegetation Index (NDVI), which was derived from the atmospherically corrected bi-directional surface reflectance data obtained from NASA Land Processes Distributed Active Archive Center’s MOD13Q1 products. The data, which had been masked to exclude water, clouds, heavy aerosols, and cloud shadows, was used to compute NDVI between 9 May and 30 September. The NDVI was calculated using the following equation.

where ρNIR is the reflectance in the near-infrared band (841–876 nm), ρRed is the reflectance in the red band (620–670 nm).

Soil DNA was extracted from 0.5 g of frozen soil samples using the FastDNA SPIN Kit for Soil by MP Biomedicals. The quantity of the extracted DNA was determined using a NanoDrop 2000 Spectrophotometer from Bio-Rad Laboratories. Polymerase chain reaction (PCR) was performed using specific primer sets for prokaryotes (515F and 806R) and fungi (ITS1F and ITS2) targeting the ITS1 region. The resulting amplicon was sequenced using the Illumina MiSeq platform and generated approximately 250 bp paired-end reads. The raw sequencing data was analyzed using the QIIME2 software (version 2020.6; Bolyen et al., 2019). The DADA2 plugin in QIIME2 (Callahan et al., 2016) was used to filter low-quality sequences and eliminate chimeras, producing amplicon sequence variants (ASVs). These ASVs were classified using the QIIME2 naive Bayes classifier (Bokulich et al., 2018), which was trained on 99% operational taxonomic units from the SILVA rRNA database (v 132) (Quast et al., 2012) for prokaryotes and the UNITE database for fungi (Nilsson et al., 2019). The raw sequence results were presented in Supplementary Table S1.

To ensure a more stringent analysis, ASVs with small counts in a limited number of samples were removed given they are likely the results of sequencing errors or low-level contamination. The analysis was based on filtered data, which only included ASVs that occurred in at least 20% of the total samples and had at least 4 counts per sample. The diversity of both the prokaryotic and fungal communities was calculated using the Shannon index for each individual sample. The community dissimilarity, or beta diversity, was estimated using the Bray-Curtis index and was visualized using principal coordinate analysis (PCoA). The statistical differences were calculated using a permutation analysis of variance (PERMANOVA). All of these analyses were performed using R 4.0.1 with the vegan package (Oksanen et al., 2013).

The co-occurrence networks in previous studies were often constructed using subjective thresholds, leading to a lack of objectivity. To address this issue, a random matrix theory (RMT)-based approach was employed to construct prokaryotic and fungal co-occurrence networks objectively using the Molecular Ecological Network Analyses Pipeline¹ and the constructed networks were visualized using Gephi². The network analysis followed the protocols established in studies by Deng et al. (2012) and Zhou et al. (2010). The network’s topological features were evaluated to determine its complexity. Keystone nodes were identified based on hub nodes including network hubs, connectors, and module hubs, which were classified based on within-module connectivity (Zi) and among-module connectivity (Pi) (Olesen et al., 2007).

Cohesion, a metric based on abundance-weighted pairwise correlations, was calculated using the protocol outlined by Herren and Mcmahon (2017). The metric is calculated as follows:

where m is the total number of taxa in a community and i is the taxa in each sample.

This protocol measures the interconnectedness of a community using a Pearson correlation matrix. The strength of each pairwise correlation was verified using a null model. The average positive and negative null model-corrected correlations were calculated for each sample to obtain a connectedness matrix. Finally, positive and negative cohesions were calculated for each sample using the above formula. Cohesion reflects the degree of cooperative behavior or competition within a community and can serve as a proxy for the strength of biotic interactions if taxa are subject to similar environmental drivers and influenced differently by species interactions. The community cohesions were then used for analyses to determine βNTI and stability measurements.

We adopted a null-model-based approach to assess the role of niche versus neutral processes in shaping microbial communities. This approach has been widely used in previous studies (Stegen et al., 2015; Zhou and Ning, 2017). The Mantel correlogram (using the “mantel.correlog” function in the R package “vegan”) was employed to estimate the correlation between niche differences (calculated as the pairwise Euclidean distances of the environmental optima of ASVs) and phylogenetic distances. The stronger signal observed at short phylogenetic distances (as shown in Supplementary Figure S2) indicated the potential for quantifying phylogenetic turnover using βNTI.

The standardized effect size of this correlation, known as the beta nearest taxon index (βNTI), was calculated as the difference between the observed value and the mean of the null distribution, normalized by its standard deviation (Stegen et al., 2013, 2015).To further distinguish the processes driving community assembly, we calculated the Bray–Curtis based Raup-Crick metric (RCbray). This metric is based on the difference between the observed Bray–Curtis dissimilarity and its null distribution (Dini-Andreote et al., 2015; Stegen et al., 2015). Values of |βNTI| <2 and RCbray > 0.95 indicate dispersal limitation, |βNTI| <2 and RCbray < −0.95 indicate homogenizing dispersal, and |βNTI| <2 and |RCbray| < 0.95 indicate drift. We estimated the relative importance of each process by calculating the fraction of values of βNTI and RCbray within different categories.

Spatial variables are determined by distance-based Moran’s eigenvector maps (dbMEM), which were called principal coordinates of neighbor matrices (PCNM; Borcard and Legendre, 2002), is a standard method for partitioning the effects of the space in ecological studies. First of all, a matrix of dbMEM variables was constructed based on latitude and longitude with the adespatial package (Dray et al., 2017).

The Moran’s I test (Moran, 1950) was utilized to assess the spatial correlation, and only eigenfunctions from dbMEM that showed a positive correlation were included in further analysis (Borcard et al., 2011). Redundancy between Hellinger-transformed prokaryotic and fungal community data (ASVs) and geographic coordinates were used to determine if there are significant linear trend. If there were significant linear trends, detrended analysis was conducted because linear trends masked all other spatial structures during variation partitioning analysis (VPA). The significance of the coefficients of the detrended matrix regression on the dbMEM matrix was tested by conducting 999 permutations of the residuals (Borcard et al., 2011). The selection of significant dbMEM variables was performed through a forward selection process, based on the 999 Monte Carlo permutation procedure of the residuals, using the “packfor” package (Blanchet et al., 2008). The variation partitioning analysis was then conducted using the significant dbMEM variables, in order to determine the relative impact of dispersal limitation on the dissimilarity of microbial communities.

The driving factors of microbial community dissimilarity can be divided into three categories, environmental variables, spatial variables (significant dbMEM) and biotic interactions (network topological features). The environmental characteristics included soil properties (i.e., pH and SWC), plant index (i.e., NDVI) and climatic factors (i.e., MAT and MAP). Random forest analysis was conducted to estimate the explained variations in driving factors to prokaryotic and fungal alpha diversity across two biomes based on R² with the randomforest package (Liaw and Wiener, 2002). The significance and impact of each predictor were determined based on the increase in mean square error (MSE) using 999 permutation tests with the rfpermute package (Archer, 2016). Meanwhile, the best regression model based on AIC was used to select the primary factors with the stats package. Then, the line regression model was used to determine their relative influence on microbial Shannon diversity. The distance-decay relationship was measured using ordinary least-squares regressions and Mantel tests between geographic distances and community dissimilarities, and the geographic distances among sites were calculated with the geosphere package (Hijmans et al., 2017). The connection between the variability in the environment and the differences in the microbial communities was calculated using a linear model analysis based on distances with the ecodist package (Goslee and Urban, 2007). The significance of the slope in the relationship between the distance and dissimilarity was evaluated through 999 permutations of residuals. Finally, the relative impact of environmental factors, spatial variables, and biological interactions were determined through a variation partitioning analysis (VPA), which was performed using the vegan package (Oksanen et al., 2013).

As shown in Table 1, most of the 26 environmental variables significantly differed across the three regions and two habitats. JQ had the lowest MAP and MAP-G but the highest PET, PET-G, AI and AI-G. HX had the lowest MAT and MAT-G. No significant differences in NDVI were observed among three regions. The altitude of HX was significantly higher than that of the other regions. HX had the highest soil nutrient contents including SWC, SOM, TN, NH₄⁺, AP and AK, while ET had the highest NO₃⁻ contents. Soil textures also differed among three regions. ET had the largest percentage of sand (80.2 ± 8.38%) and silt (16.45 ± 6.97%) while the lowest percentage of clay (3.35 ± 1.53%). The highest EC, salinity and pH were observed in HX.

Artemisia and Poaceae had significant differences in MAP and MAP-G while no differences in MAT, MAT-G, PET-G, PET, AI-G and AI. The NDVI was significantly higher in Poaceae (0.5 ± 0.1) than in Artemisia (0.22 ± 0.08). Poaceae had higher NDVI, SWC, EC, Clay, Silt, SOM, NH₄⁺, and AK contents. Artemisia had higher KS, sand, TN, NO₃-and AP contents. No significant differences in pH were observed between habitats.

The most abundant prokaryotic phyla were Actinobacteria, Proteobacteria, Acidobacteria and Crenarchaeota, which showed distinct distributions among regions but showed no significant differences between biomes (Supplementary Figure S1B). Acidobacteria and Crenarchaeota were more abundant in HX, while the relative abundance of Gemmatimonadetes was significantly higher in ET (Supplementary Figure S1A). Fungal communities also exhibited compositional differences in regions but no differences in habitats. JQ harboured a relatively higher abundance of Ascomycota but had the least abundance of Mortierellomycota, which was highest in ET.

Shannon’s diversity measurements indicated that prokaryotic communities in the ET were the highest diverse, followed by JQ and HX (Figure 2A). Poaceae had a greater α-diversity than Artemisia (Figure 2C). We found similar patterns in the fungal matrix, and ET had a greater α-diversity in comparison with other regions (Figure 2B). Shannon’s diversity in Poaceae was greater than that in Artemisia (Figure 2D).

PCoA results showed that prokaryotic communities in HX were far away from those in the other samples (Figure 2E), while no distinct separation among the three fungal communities was observed (Figure 2F). We observed no apparent separation of prokaryotic and fungal communities between Poaceae and Artemisia (Figures 2G,H). The results of the PERMANOVA test showed that the composition of both the soil prokaryotic and fungal communities were significantly different across the three regions (R² = 0.15, p < 0.001 for prokaryotic communities; R² = 0.06, p < 0.001 for fungal communities) and between two habitats (R² = 0.02, p < 0.001 for prokaryotic communities; R² = 0.06, p < 0.001 for fungal communities).

The null model revealed distinct assembly patterns between types of microorganisms but similar patterns between regions or habitats. As shown in Figure 3, prokaryotic community assembly in all regions and habitats was mainly regulated by homogeneous selection, with some impact from ecological drift. For fungi, community assembly was largely influenced by ecological drift. This was indicated by the βNTI and Rcbray values, with most of the βNTI values for prokaryotes less than−2 and most for fungi between−2 and 2, and the absolute values of Rcbray less than 0.95 for both prokaryotes and fungi (Figure 3). These findings suggest that deterministic processes were the primary drivers of prokaryotic community assembly, while stochastic processes played a more prominent role in fungal community assembly.

We found distinct differences in network structures between prokaryotic and fungal communities (Figure 4A), whereas relatively few differences were found among regions and between biomes (Supplementary Table S2). Meanwhile, we observed more edges and higher average degrees in prokaryotic networks while more modules in fungal networks, indicating that prokaryotic networks were more connected and fungal networks were more modularized (Supplementary Table S2). We then explored the biotic interactions of soil microorganisms by establishing cross-kingdom co-occurrence networks consisting of prokaryotes and fungal taxa based on correlation matrix. These results showed that fungi-prokaryotes network in Artemisia was more connected and clustered than that in Poaceae.

We found that most network attributes played crucial roles in both prokaryotic and fungal α-diversity, and network vertices had the largest observed effect (Figure 5A). Of the environmental factors, only salinity and MAT affected both prokaryotic and fungal α-diversity. Clay and EC had significant effects on the prokaryotic diversity of Artemisia and Poaceae, respectively. SOM and KS influenced the fungal diversity of Poaceae and Artemisia, respectively (Figure 5A). The species diversity between prokaryotes and fungi had no significant relationship, whereas prokaryotic network vertex increased linearly with fungal network vertex in Artemisia (Figure 5B). Based on the results of the best regression model with the lowest AIC, we found that network attributes including vertex, positive cohesion, negative cohesion were the best predictors for microbial Shannon diversity (Figure 5C).

The results of Mantel Tests showed biotic interactions (Inter) had higher correlations than environmental matrix (Env) and geographic distance (Disp) with prokaryotic community dissimilarity across all three regions and two habitats (Figure 6A). In the entire dataset and in the Artemisia habitat, the dissimilarity of the fungal community showed stronger associations with geographic distance, whereas the strongest correlations were observed between biotic interactions and fungal community dissimilarity across three regions and in the Proaceae habitat (Figure 6A). Furthermore, the results of the Variance Decomposition Analysis showed that biotic interactions played a greater role in shaping both the prokaryotic and fungal communities than environmental factors and geographic distance, both in the overall dataset, habitats of Poaceae and Artemisia (Figure 6B) and three regions (Supplementary Figure S3).

The results of Partial Mantel Test showed that biotic interactions showed the strongest correlation with both bray–curtis dissimilarity (bray) and phylogenetic turnover (βNTI) of prokaryotic community, followed by environmental factors and geographic distance. Fungal βNTI showed the highest correlation with biotic interactions but fungal bray-curits dissimilarity had higher correlations with geographic distance than biotic interactions (Table 2).

The null model analysis based on βNTI and Rcbray could distinguish the relative contributions of the deterministic process and stochastic process in governing the microbial communities (Stegen et al., 2013, 2015). In this study, most of the βNTI values for prokaryotes were less than−2 and most for fungi were between−2 and 2, suggesting that deterministic processes played a stronger role in prokaryotic community assembly, whereas fungal community assembly was structured by stochastic processes (Figure 3). |Rcbray| < 0.95 further indicated that stochastic processes were dominanted by ecological drift (random proliferation or death of microorganisms). Such differential responses across taxonomic types have also been observed in various habitats, such as soil (Powell et al., 2015), glaciers (Jiang et al., 2018), sediment (Zhao et al., 2022) and freshwater (Logares et al., 2018). Despite significant differences in soil physicochemical properties and spatial distance (Table 1), microorganisms belonging to the same groups of soil microorganisms showed similar assembly modes in different geographic regions or habitats (Figure 3). These results support our first hypothesis that taxonomic type plays a greater role than habitats types and geographic distance on soil microbial community structures.

One possible explanation for the differential response of prokaryotic and fungal communities to geographic distance and habitat types is that they have distinct sensitivities to environmental selection. In addition, increasing environmental heterogeneity can result in greater dissimilarity among microbial communities with increasing geographic scale. The VPA and Partial mantel test showed that changes in habitat environments have a greater impact on community turnover in prokaryotes compared to fungi (Figure 6; Supplementary Table S3), suggesting that prokaryotes are more susceptible to environmental changes than fungi. Moreover, the differences in dispersal capacity of the two groups of microorganisms could also foster the distinct assembly models. Propagule sizes of fungi are typically within the range 5–50 μm diametre (Ingold, 1971) while those of bacteria are usually within 0.2–20 μm diametre (Young, 2006), organisms with smaller propagule size and larger number are typically thought to be easier to spread. Mantel test between distance and community structure revealed the differences of Prokaryotes and fungi in dispersal capacity (Figure 6A). However, large spatial distance could weaken the difference in dispersal capacity in soil between prokaryotes and fungi. A survey reported the limited impact of dispersal capacity in determining the distributions of soil bacteria at global scale (Bahram et al., 2018). The weaker correlations between geographic distances and both prokaryotic and fungal dissimilarities further supported less effects of dispersal limitation on soil microbial community assembly at region or bigger scale (Figure 6A).

Ecological drift, a stochastic process referring to random births and deaths in a population (Vellend, 2010), is expected to be more significant with decreasing population size because random demographic events play a significant role in smaller populations (Fodelianakis et al., 2021). Fungi with a smaller population size relative to prokaryotes hence tend to be affected by the greater role of ecological drift. The results based on null model analysis also suggested that fungal community assembly was dominated by ecological drift. The large part of the unexplained variance of fungal community turnover in VPA may also be due to the unmeasurable stochastic process.

A recent global-scale study reported that environmental filtering was the predominant driver of soil prokaryotic community turnover (Fierer and Jackson, 2006). We also observed that the community turnover of prokaryotes showed a more significant relationship with the aggregate environmental matrix than that of fungi (Figure 5B). VPA revealed that the environmental matrix explained more variance in community dissimilarity of prokaryotes than that of fungi, further demonstrating the great role of environmental selection in prokaryotic community assembly (Figure 6A). These results further indicate that deterministic processes (environmental filtering) have more impacts on soil prokaryotic community, while soil fungal community were mainly structured by stochastic processes.

The results of the Mantel analysis indicated that salinity was the strongest factor affecting both the βNTI and β-diversity (Supplementary Table S3), which is consistent with results found in desert (Zhang et al., 2019) and lake ecosystems (Yang et al., 2019). The OLS model indicated that higher salinity levels were associated with a reduction in species diversity. This negative relationship possibly due to the heightened extracellular osmolarity resulting from excessive salt concentrations (Oren, 2011; Rath and Rousk, 2015), which could cause the death or inactivity of taxa that cannot tolerate osmotic stress (Pontarp et al., 2012). Furthermore, significant correlations were found between the prokaryotic community and multiple environmental variables (e.g., MAT and AP, Figure 5A; Supplementary Table S3), suggesting that the formation of prokaryotic communities is influenced by various niches. However, it should be noted that a considerable amount of variation in the soil microbial community observed in our study was not accounted for. Possible reasons for this include unmeasured environmental variables, the limited scope of our sampling, and potential biotic interactions (Jiao et al., 2020).

Notably, evaluations of the relative contributions of determinism to community assembly largely concentrate on the set of environmental factors. However, the dominant processes refer not only to environmental selection but also to all ecological forces driving community turnover, such as biotic interactions (Singh et al., 2009). The results of the null model showed that homogeneous selection was dominant in the prokaryotic community assembly, which demonstrated that the biotic interactions existed (Danczak et al., 2018). However, elucidating microbial spatial–temporal distribution from the perspective of interspecific interactions is a tremendous challenge because it is more difficult to directly quantify the interactive patterns of microorganisms relative to macroorganisms. A typical way of including species interactions into an explanation matrix is to apply correlation-based cooccurrence network analysis (Layeghifard et al., 2017).

Random forest analysis revealed that network topological parameters had better prediction performance in both prokaryotic and fungal alpha diversity than environmental variables (Figure 5A). Linear regression analysis found that the number of network vertices and positive cohesion could explain more variance in microbial Shannon diversity (Figure 5C). These findings indicated that network topological features could be an effective proxy for biotic interactions and that biotic interactions play a critical role in driving microbial alpha diversity.

Biotic interactions have been reported a predominant factor in the β-diversity of diazotrophic and bacterial communities in paddy soil (Gao et al., 2019). In this study, the network topological feature matrix had stronger correlations with prokaryotic and fungal β diversity (Figure 6A), which indicates biotic interactions may also be the primary driving in structuring soil prokaryotic and fungal community composition. In addition, we found network attributes had stronger correlations (Figure 6A; Table 2) and explained more variance of both prokaryotic and fungal community structure dissimilarity and phylogenetic turnover than environmental factors and geographic distance (Figure 6B). These results suggested that taxa-taxa interactions may play a more important roles than environmental filtering or dispersal limitation in both soil prokaryotic and fungal community assembly processes.

Community complexity was investigated by relating a recently published measure of cohesion of biotic interactive network to ecological structuring processes (Herren and Mcmahon, 2017). Negative cohesion of microbial co-occurrence network has been reported to be significantly correlated with the β-diversity of bacteria (Herren and Mcmahon, 2017). A study in aquifers found that microbial communities with more negative cohesion values experienced lower turnover and were more likely to regulated by homogenizing selection, whereas less complexity communities experienced higher turnover and susceptibility to stochastic processes (Danczak et al., 2018). In this study, we also found significantly positive correlations between negative cohesion and βNTI values of prokaryotic communities in arid ecosystem (Supplementary Figure S4), indicating not only the connectivity but also the complexity of biotic interactions could affect the microbial community turnover and community assembly processes.

A previous study reported that fungal richness could affect the relative contributions of deterministic and stochastic to bacterial community assembly processes (Jiao et al., 2021). Complex interactions between soil fungi and prokaryotes in the cross-kingdom network were observed in this study (Figure 4). These results may indicate that the assembly mechanism of prokaryotic and fungal communities was influenced not only by intra-kingdom interaction but also by the biotic interactions across taxonomic types (Duan et al., 2022). Despite no significant correlation between the species diversity of fungi and prokaryotes in Artemisia habitat, we found significant positive correlations between their vertex (Figure 5B), which implies the potential facilitation or symbiosis between prokaryotic and fungal taxa (Duan et al., 2022). Therefore, estimation based only on intra-kingdom biotic interactions may introduce bias, and incorporation of cross-kingdom interaction could facilitate the understanding of community assembly, especially in networks with intense inter-kingdom interactions.

It is noteworthy that the effects of biotic interactions were estimated based on the microbial co-occurrence network, which is only a putative species interaction network yielding statistical associations between taxa (Carr et al., 2019). It cannot be too prudent to infer the actual microbial interactions and their ecological meaning. Despite unavoidable defects, the spatiotemporal dynamics of co-occurrence networks hold the potential to affect community assembly and topological features could serve as an effective proxy to estimate the relationships between biotic interactions and ecological processes.

Together, our study revealed the potential of microbial biotic interactions to serve as predictors or interpreters for prokaryotic and fungal α-and β-diversity, built the linkage between diverse microbial co-occurrence network topological features and mechanisms underlying microbial community assembly, and highlighted the importance of biotic interactions including intra-and cross–kingdom interactions, in regulating microbial diversity and microbial community assembly processes across terrestrial ecosystems.