A comprehensive sampling of biocrusts was collected across seven major deserts and 3,500 km in northern China. A total of 200 samples were collected from 13 plots, including 140 A, 24 C, and 36 M-dominated biocrusts. Each sample plot was considered an SP, as shown in Figure 1. Each sample plot had received no rainfall in the past 72 h and kept 0.2 m away from the shrubs. The biocrust sample covered by a single biotype was selected for each successional stage. The biocrusts and attached subsoil were gathered with a shovel and preserved in sterilized plastic Petri dishes to ensure the integrality and then transported to the laboratory within 12 h. To ensure the non-redundancy and representativeness of the sequencing results, the macroscopic moss plants (but not their protonemata) were removed from M biocrusts (Li and Hu, 2021).

Total genomic DNA was extracted from biocrust samples using PowerSoil^® DNA Isolation Kit (MOBIO, United States). 16S rDNA (524F10extF/Arch958RmodR, V4–V5 region) primers were used for amplification. Sequencing was performed on the Illumina MiSeq PE300 platform (Illumina, San Diego, CA, United States). The raw sequence data were uploaded on NCBI under BioProject PRJNA730649. The acquired sequences were filtered for quality control following standard procedures (Shanghai Majorbio Bio-Pharm Technology Co., Ltd.). The operational taxonomic unit (OTU) was generated from defined representative sequences, which set clustering at 97% similarity. Ribosomal database project (RDP) classifier was used for taxonomic annotation of representative sequences based on an identity threshold of 0.7 in the SILVA 132 database for archaea (16S_archaea).

Microenvironments included water content (WC), soil texture (%), thickness (TH., mm), chlorophyll a (Chl a, μg⋅g^–1), scytonemin (Scyt., unit⋅mg^–1 FW), extracellular polysaccharide (EPS, mg⋅g^–1), bacteriochlorophyll a (BChl a, μg⋅g^–1), variable fluorescence/maximal fluorescence (Fv/Fm), pH, ORP, salinity (μmol⋅g^–1), NH₄⁺ (μmol⋅g^–1), PO₄^3– (μmol⋅g^–1), total nitrogen (TN, g⋅kg^–1), total phosphorus (TP, g⋅kg^–1), total organic carbon (TOC, g⋅kg^–1), soil alkaline protease (ALPT), soil-β-glucosidase (β-GC), and soil alkaline phosphatase (ALP). Macroclimates included mean annual precipitation (MAP), aridity index (AI), mean annual sunshine duration (MASD), mean annual temperature (MAT), altitude (Alt.), and wind speed (WS). These variables were determined by previous methods (Li and Hu, 2021).

To study the driving forces on the difference of archaeal communities varied with the number of SPs, the combined meta-community dataset containing the different number of sample plots was derived by the method of combination. First, there were 13 sample plots, and each was considered an SP. We calculated C133, C134, C135, C136, C137, C138, C139, and C1310, corresponding to three to ten SPs. Each combined dataset included at least one of the Fukang, Zhangye, Minqin, and Shapotou sample plots and all three successional biocrusts (Figure 1). 202, 589, 1161, 1632, 1680, 1278, 714 and 286 combined communities were generated under the three to ten SPs respectively. This step enabled us to measure the influence of successional factors on the variation of combined meta-communities. Second, the permutational multivariate analysis of variance (PERMANOVA, permutation = 999) was used to calculate the explanations (R²) of successional stages (i.e., successional factor), geographical locations (i.e., spatial factor), and each environmental variable (i.e., environmental factor) for each combined meta-community variation-based Bray–Curtis and UniFrac-Weighted dissimilarity. Third, the geographic distance of a combined meta-community dataset was defined by the distance between the farthest sample plots in the given dataset. Subsequently, the linear regression was conducted between the explanation of these three factors and geographic distance (*<0.05, ^**<0.01), and the slope was used to indicate the driving forces of the three factors affecting the biogeographic patterns. Note that the explanation of environmental factors for the given combined meta-community was represented by a weighted mean explanation of all environmental factors to avoid Simpson’s Paradox error. The phylogenetic tree was constructed by IQ-TREE (version 1.6.8¹) based on the maximum likelihood method with default procedures and was used to calculate the UniFrac-Weighted dissimilarity.

To analyze the characteristics of archaeal community structure and its covariant relationship with environments, the redundancy analysis (RDA) was performed at a family level. Then, the forward selection procedure was used to screen out significant environmental variables (Blanchet et al., 2008). The analysis of similarities (ANOSIM) (permutation test = 999) was used to examine community differences among successional stages and between subgroups. Linear regression analysis (*p < 0.05, ^**p < 0.01) was used to assess the goodness of fit (R²) for the first axis coordinates of non-metric multidimensional scaling (NMDS) with environmental factors (standardization) in each successional stage. To explore the correlation between community variation and environmental variables, a mantel test analysis was performed with community matrices at OTU level and environmental factors. ANOVA was used to analyze the significant differences of α and β diversity between three successional stages [least significant difference (LSD), *p < 0.05, ^**p < 0.01]. All of the above results were visualized by box and whisker (10–90%). We referred to Stegen’s method (Stegen et al., 2013) to quantify the relative importance of ecological processes. β-Nearest-taxon-index (βNTI) values of >2 or <−2 indicated heterogeneous selection or homogeneous selection, respectively, whereas if the absolute values of βNTI were <2, stochastic processes would play an important role. Then, the Raup–Crick metric using Bray–Curtis dissimilarities (RC_bray) values <−0.95 pointed to a community assembly governed by homogeneous dispersal. In contrast, dispersal limitation could generate RC_bray values >+0.95. Drift is only expected when RC_bray values were between −0.95 and +0.95. We inferred a co-occurrence network for the archaeal community at the OTU level using the default settings on the MENA website² and generated characteristics of topological structure. We identified key and peripheral species in the network by the within-module connectivity (Zi) and among-module connectivity (Pi) based on Deng’s method (Deng et al., 2012).

The R environment (version 3.6.2³) was used for relevant statistical analyses.

The α diversity of archaeal communities in A and M crusts was significantly higher than that in C crusts (Figure 2A). The results of ANOSIM showed that there was no significant difference in the archaeal community among the three successional stages (Figure 2B). In contrast, in the same successional stage, the differences of archaeal communities in A and M crusts were larger than that in C crusts and the difference in A crusts has exceeded the difference between succession. The results of ecological processes showed that the archaeal communities were mainly assembled by drift in A and M crusts, while they were mainly dominated by homogenous dispersal in C crust (Figure 2C). The co-occurrence patterns demonstrated that mutual exclusion was dominant in A and M crusts, while coexistence was dominant in C crust (Figure 2D). In summary, the assemblages mainly governed by drift had higher α diversity, community differences, and mutual exclusion in A and M crusts than in C crusts.

The biogeographic patterns based on both taxonomic and phylogenetic variations, which are explained by successional, spatial, and environmental factors, were examined (Figure 3 and Supplementary Figure 1). In terms of slope direction, the explanation of spatial and environmental factors increased with distance, but successional factors decreased. At the taxonomic variation facet, the order of slope absolute value was spatial > successional > environmental, and the slope absolute value of spatial and successional factors decreased with the number of SPs. In addition, based on phylogenetic variation, the order was spatial > successional > environmental at three and four SPs. In contrast, with five to ten SPs, the order was spatial > environmental > successional, and the slope absolute value of spatial and successional order decreased. In summary, spatial and environmental factors were positive driving forces at both taxonomic and phylogenetic facets of archaeal community variation, while successional factors were negative in three to ten SPs. The impact of the spatial factor was the most significant, but it differed at both taxonomic and phylogenetic variation facets.

Analysis of similarities and RDA were, respectively, used to examine differences between the archaeal community structure and significant environments at the family level. The total variation of the archaeal community was explained most by AI, followed by salinity and PO₄³⁺ (Supplementary Table 2). Two hundred samples were clustered into two subgroups: subgroup1 (SG1) and subgroup2 (SG2) (Figure 4A; ANOSIM: R² = 0.98, p = 0.001). Meanwhile, the correlation of community structure with environments between SG1 and SG2 sharply contrasts (Figure 4A). Specifically, the community structure in SG2 had a negative correlation with Fv/Fm, NH₄⁺, oxidation-reduction potential (ORP, mv), MAP (mm), Alt. (m), and WC (%), while had a positive correlation with PO₄³⁺, salinity, silt, WS (m⋅s^–1), pH, AI, MAT (°C), and MASD (h). These correlations are reversed in SG1. The results of the community structure at the family level indicated that Nitrososphaeraceae in Thaumarchaeota was dominant in SG1 (Figure 4B), followed by unclassified_d_archaea. In our study, Halococcaceae and Haloferacaceae were mainly Haloarchaea, which dominated in SG2. In summary, the abundance of various taxa varied across 200 archaeal communities. Nitrososphaeraceae dominated the archaeal communities in SG1, and their distribution pattern was positively associated with non-salinity-related environmental variables. In contrast, Haloarchaea dominated SG2, and the distribution was positively associated with salinity-related environmental variables.

We, respectively, quantified the relative roles of ecological processes in SG1 and SG2 and discovered that SG1 assembly was primarily influenced by drift alone (Figure 5), followed by a homogenous selection, while SG2 assembly was dominated by the homogeneous selection, which was significantly larger than drift. Comparatively speaking, the deterministic process had a greater impact on SG2’s community assembly than on SG1.

We derived co-occurrence networks for SG1 and SG2 (Figure 6 and Table 1). In terms of topological structure, mutual exclusion dominated the relationships of both subgroups but was stronger in SG1. However, SG2 contains more nodes, a larger number of edges per node, an average degree, and a longer average path distance than SG1. In terms of co-occurrence species, unclassified_d_archaea and Nitrososphaeraceae dominated in SG1, while Haloarchaea was barely involved. However, unclassified_d_archaea, Nitrososphaeraceae, and Haloarchaea all dominated in SG2. The relative proportion of peripheral taxa was slightly higher than that of key taxa in SG1 (55.71 vs. 44.29%), whereas the contrast was more evident in SG2 (92.46 vs. 7.54%). In summary, archaea co-occurred in a mutually exclusive pattern. The SG1 network had a higher proportion of competitions and key taxa than the SG2 network, although the network of SG2 was more complex.

Many ecological effects have been suggested to understand biodiversity changes including “mass effects” and “drift,” which vary widely in importance depending on the environments (Leibold et al., 2004; Li and Hu, 2021). In particular, the effects of habitats on community heterogeneity were more emphasized in archaeal communities (Jiao et al., 2019). Our results showed that there were more macroclimate factors than microenvironments in factors that had greater explanation compared to succession factor (Supplementary Figure 1). This may result in a consequence of that the archaeal community variations in the same successional stage exceed the variations between successional stages (Figure 2B), which was quite distinct from the successional convergence of prokaryotic community in biocrusts (Xu et al., 2020). In other words, the criterion of classification for different successional stages of biocrusts based on the dominant cryptogamous plants was not applicable to the archaeal communities’ development. One possible explanation for this fact is that the archaea with low abundance in biocrusts were more vulnerable to the effects of drift (Pedros-Alio, 2006). On the other hand, archaea with unique niche were facilitated to exchange among unknown patchy habitats in the way of dispersal under the flat and windy geomorphic characteristics of northern China (Wilkinson et al., 2012). Further combining the results of biodiversity, assembly changes, and co-occurrence patterns (Figure 2), it can be inferred that archaeal communities were assembled by heterogeneous taxa driven by drift in A and M crusts, resulting in a higher α, β diversity, and mutual exclusion. In contrast, the archaeal communities in C crusts were assembled by homogeneous taxa driven by “mass effects,” which showed a low α, β diversity, and coexistence. In addition, the distributions of archaea in A and M crusts were mainly affected by microenvironments, and that in C crust were mainly controlled by macroclimates (Supplementary Table 1). This suggested that microenvironment and macroclimate variables may control the ecological processes of archaeal community assembly (Dini-Andreote et al., 2015; Zhang et al., 2016). In summary, based on the results of biodiversity, assembly, and co-occurrence patterns, we proposed that it was more suitable to classify archaeal communities by dominant taxa rather than by succession for understanding the ecological laws. Meanwhile, our results also suggested that patchy habitats played an important role in archaeal community variation.

It is crucial to determine the biogeographic patterns and their driving forces in microbial ecology. Considering that microbial community heterogeneity ranges from a few centimeters (Chu et al., 2016) to thousands of kilometers (Crowther et al., 2019), it implied that the biogeographic patterns may depend on the number of SPs considered. Commonly, the slope of factors explaining community heterogeneity with distance was used to characterize the driving force. In addition to spatial factors (Ramette and Tiedje, 2007; Martiny et al., 2011; Wang et al., 2017), it was used to examine environmental filtering effects (Chen et al., 2019; Jiao et al., 2020) and the regulation with succession (Gotelli et al., 2017). In our study, the driving force of spatial factors was greater than that of environmental and successional factors (Figure 3), and it also had the greatest explanation for archaeal community difference (Supplementary Figure 1) under the number of three to ten SPs, which could be because the number of isolated SPs of archaeal communities was higher than expected. The pronounced effects of spatial factors include not only geomorphic but also the comprehensive effects of unmeasured transient environments (Ouyang and Hu, 2017), niche conservation, dispersal rates, and historical events (Hendershot et al., 2017). The spatial effects are amplified at large scales (Kramer-Schadt et al., 2013; Li and Hu, 2021). However, another study reported that archaeal community heterogeneity was primarily driven by environments in the agricultural ecosystem (Jiao et al., 2019). The difference in driving factors could be due to the extinction of niche conservative archaea under artificial periodic interference in the agricultural ecosystem. This situation was rarely observed in dryland (Figure 3), although the strength of environmental filtering was always maintained at a high level. However, the fluctuation of macroclimates was more usual and irregular compared to the microenvironments, as indicated by the result that AI explained 48% of the total community variation (Supplementary Table 2), which sparked our interest in comparing macroclimate and microenvironment variables. Our results revealed that macroclimates explained more variation in archaeal communities than microenvironments (Supplementary Figure 2). In recent studies, succession as a factor that converged community differences has been revealed in bacterial and eukaryotic communities in biocrusts (Xu et al., 2020; Li and Hu, 2021), and comparable but lesser effects of succession have been observed in archaeal communities (Figure 3). Therefore, the ability of various microorganisms to self-organize coincided with temporal trends in the biocrust ecosystem. In summary, spatial factors were critical in improving the archaeal community heterogeneity, and succession was the only factor that converged community differences.

Microbial taxa shared among habitats would be limited by their dispersal abilities and spatial distance (Liao et al., 2017; Mo et al., 2018; Li et al., 2021), and less may be shared in our sampling transect with obvious altitude gradients along the meridional direction. Therefore, we assumed that the number of isolated SPs considered could alter biogeographic patterns (Lennon and Jones, 2011; Kramer-Schadt et al., 2013). According to our results (Figure 3), successional and spatial factors drove the most apparent biogeographic patterns under the number of three SPs at the taxonomic facet, highlighting the considerable heterogeneity of archaeal assemblage composition suffering from microhabitat changes (Jiao et al., 2019, 2021). In contrast, at the phylogenetic facet, five SPs were the turning point where the slope of spatial factors began to decrease (Figure 3), and it also indicated the tipping point where the influences of community homogeneity began to emerge, which was mostly contributed by the widely distributed Nitrososphaeraceae (Figure 4A). Meanwhile, increasing homogeneity suppressed the successional driving force (Xu et al., 2020), explaining our results that weaker biogeographic patterns driven by successional factors were exhibited with an increasing number of SPs considered. Overall, the spatial driving force for taxonomic and phylogenetic variation was not synchronized with the increasing SP, but the environmental and succession driving forces were synchronized. In addition, the increasing community homogeneity from the widespread Nitrososphaeraceae suppressed the successional driving force.

Most archaea inhabit marine sediments and 8–5,000 m deep underground, accounting for approximately 20% of the total microbial community, while only 1% exists in the dryland topsoil with frequent WS changes, high light intensity, and temperature (Bates et al., 2011; Bar-On et al., 2018). According to our results, archaeal communities were dominated by Nitrososphaeraceae and Haloarchaea, whose distribution patterns were affected by non-salinity-related and salinity-related variables, respectively. Nitrososphaeraceae abundance demonstrated a negative correlation with MAT and MASD (Figure 4A), which was consistent with their non-thermophilic attribute (Tourna et al., 2011). It was also understandable that the distribution of SG1 can be explained by NH₄ to a certain extent (Figure 4A), because NH₄⁺ was inclined to exist as NH₃ to participate in the ammonia oxidation process under leaning alkalinity environment of biocrusts (Yuan et al., 2019). Considering the obvious altitude gradients along the meridian of China, Alt. should contain more complicated meanings besides unmeasured moisture. In addition, MAP and WC also altered Cyanobacteria composition in dryland topsoil, affecting the organic carbon quality (Yuan et al., 2016). When this was combined with our result that Nitrososphaeraceae abundance exhibited a positive correlation with MAP and WC (Figure 4A), we concluded that Nitrososphaeraceae was associated with organic carbon quality, which implied that they could be heterotrophs. Cyanobacteria, which are also major nitrogen-fixing microorganisms in biocrusts, contribute most to Fv/Fm, a characteristic of photosynthetic capacity. The distribution of SG1 was affected by Fv/Fm (Figure 4A) implying that the ammonia oxidation process could be coupled with nitrogen fixation.

Haloarchaea is a branch of non-methanogenic archaea that thrives in salinity-dependent environments, due to a “salt-in” osmoprotection mechanism (Becker et al., 2014). This can also be proved by our results of the positive correlation between their abundance and WS, MAT, MASD, and AI (Figure 4A). Moreover, they have an adaptive evolution from anaerobic to aerobic habitats (Sorokin et al., 2017). Our results indicated that the communities of SG2 distributed along with the content of silt (Figure 4), which is easy to form anaerobic microchamber with low light (Chau et al., 2011), implying a relatively older phylogenetic position. In addition, phosphorus is not only a limiting variable of terrestrial ecosystem productivity (Cleveland et al., 2013) but also a key variable in regulating the assembly of bacterial and eukaryotic communities in biocrusts (Li and Hu, 2021). The effect of PO₄³⁺ on the distribution of SG2 communities in this study implied that it acted as a valve in the process of community assembly and productivity allocation. The results of the mantel test also showed that there were significant relationships between SG2 and ALPT, ALP, and β-GC (Supplementary Figure 3), which suggested that Haloarchaea in biocrusts contributed to the soil carbon cycle by degrading extracellular proteins, carbohydrates, and straight-chain lipids (Iverson et al., 2012; Li et al., 2015).

In microbial ecology, the assembly mechanism shaped biogeographic patterns. However, much remains unknown about the archaeal community assembly. To the best of our knowledge, this is the first study to reveal archaeal assembly mechanisms in dryland topsoil habitats. Our results indicated that the dominant assembly processes in SG1 and SG2 were completely different (Figure 5). First, the dynamic community assembly processes were generally regulated by environments with the greatest heterogeneity (Logares et al., 2013; Liu et al., 2020). Considering this, WC and salinity could be considered key environments for regulating assembly in SG1 and SG2, respectively (Figure 4A). Second, the spatial scale was also a non-negligible factor in altering different microorganism assembly processes. Our result indicated that stochastic processes dominated in SG1 with larger spatial scale, whereas deterministic processes dominated in SG2 with smaller spatial scale (Figure 5). This situation was consistent with the bacterial community assembly (Chen et al., 2020; Liu et al., 2020; Logares et al., 2020) but not with eukaryotic communities (Chen et al., 2017; Zhang W. et al., 2018). Third, the ability of ATPase operons to horizontally transfer expands the habitat range of Nitrososphaeraceae (Wang B. et al., 2019), which may be a more biased factor to easily detect the influence of drift based on the assembly quantitative method of phylogenetic signals. In contrast, Haloarchaea has the ability of light utilization through rhodopsin (Kamo et al., 2006), which could be restricted by the microchamber with low porosity in biocrusts considered as a deterministic factor affecting assembly. Fourth, under the same size of SP (Supplementary Figure 4), the assembly patterns of meta-communities dominated by Nitrososphaeraceae were governed by drift, while that dominated by Haloarchaea were governed by homogeneous selection, which implied that the assembly patterns of SG1 and SG2 were determined by the dominant taxa rather than the size of SPs. Collectively, the comprehensive influence of environmental, spatial, methodology, and survival strategies resulted in the different assemblies in SG1 and SG2.

Many studies concerning various habitats and microorganisms revealed intricate co-occurrence relationships (Liu et al., 2020; Li et al., 2021). In biocrusts, mutual exclusion dominated in archaeal communities (Figure 6), which was completely different from the high coexistence of bacterial and eukaryotic communities (Li and Hu, 2021), but was similar to archaea in river sediments resulting from the narrowing of niche breadth due to environmental filtering (Chen et al., 2020). This could be a possible explanation for the distribution patterns in our result that are related to different attributes and environmental variables (Figure 4A). Moreover, the higher proportion of homogenous microorganisms in the SG1 and SG2 (Figure 4B) resulted in microbes occupying similar resources, let alone the Nitrososphaeraceae predominate ammonia-oxidizing process in the nitrogen cycle, which enhanced the competition by simplex ecological function. Furthermore, the high proportion of key taxa in SG1 (Table 1) could contribute to the drift in the community assembly because the competition made adaptation to a new habitat more difficult. In comparison, the competition in the SG2 network was slightly weaker with less proportion of key taxa, enabling the appearance of homogeneous selection in the community assembly. In addition, the SG2 co-occurrence network had a high proportion of Haloarchaea (Figure 6), which was related to their ability to use rhodopsin for photoautotrophy (Kamo et al., 2006), and “salt-in” strategy that provided transmembrane dynamic potential in a hypersaline environment. Furthermore, in the same size of SP, the meta-communities dominated by Nitrososphaeraceae had more proportion of mutual exclusion and key taxa than that dominated by Haloarchaea (Supplementary Table 3). In summary, we believed that the dominant taxa with different survival strategies in archaeal communities determined the co-occurrence patterns and proportion of key taxa, which was also a key to understand different archaeal community assembly patterns.