Structure, Assembly, and Co-Occurrence of Archaeal Communities in Tree Rhizosphere of the Qinghai-Tibetan Plateau

Background: Archaea are considered to be an important component of complex microbiomes, which play vital roles in mediating soil biogeochemical processes. However, little is known about the ecology of archaeal communities in the rhizosphere, especially regarding their assembly and co-occurrence patterns. Results: Here, we used high-throughput Illumina sequencing to investigate the community variations of archaea between the rhizosphere and bulk soil collected from two native alpine tree species of the Qinghai-Tibetan Plateau. Thaumarchaeota was the dominant archaeal phylum in all soils tested, while archaeal community structures in the rhizosphere signicantly differed from that in the bulk soil. Soil ammonium nitrogen, soil organic matter, available phosphorus and pH were signicantly correlated with the archaeal community structure, and the deterministic processes dominated the assembly of archaeal communities across all soils. In addition, the network structures of the archaeal community in the rhizosphere were less complex than they were in the bulk soil, and an unclassied archaeal group (Unclassied_k_norank) was identied as the keystone species in all archaeal networks. Conclusions: Collectively, our ndings suggest the structural variability, assembly processes, and co-occurrence patterns of archaeal communities in the rhizosphere of the the Qinghai-Tibetan Plateau, which further deepens our ecological understanding of the archaeal microbiome.


Background
The rhizosphere is a narrow zone of soil that tightly surrounds growing plant roots, which secrete a variable but substantial amount of photosynthesis-derived organic carbon compounds that enable the growth and metabolic activities of soil microorganisms [1,2]. Therefore, the rhizosphere has been considered to be one of the most complex interfaces in nature [3], where a variety of microorganisms drive multiple biogeochemical transformations including soil formation, carbon and nitrogen cycling [4,5]. In addition, rhizosphere microbial communities also have important effects on plant growth, health, and abiotic stress tolerance [6][7][8]. A growing number of studies have investigated the structure and assembly process of rhizosphere microbial communities, as well as their response to the selective effects of various biotic and abiotic factors [9][10][11]. However, these studies largely focused on bacteria and fungi, and little is known about the structural characteristics and driving factors of archaeal communities in the rhizosphere.
In fact, archaea have been considered a substantial component of complex microbiomes [12], and have profound interactions with bacteria, fungi and viruses in a wide range of Earth's ecosystems [13,14].
Compared to soil bacteria and fungi, archaeal communities are usually of low abundance and have less diversity [15], and they used to be thought to occur only in extreme environments [16]. Due to the rapid development of high-throughput sequencing technology, recent studies have expanded our knowledge of the biology of the archaea and have discovered their fundamental and even crucial ecological functions including methanogenesis [17], ammonia oxidation [18], hydrocarbon degradation [19], sulfate reduction [20], etc. Thus, a better knowledge of the structure, assembly and interaction of archaeal components in soil is of great importance [12].
Several studies have investigated the diversity and composition of archaeal communities in the rhizosphere. These studies are mainly limited to rice and a few wetland plants [21,22], as well as only focusing on a minority of archaea taxa such as ammonia-oxidizing archaea (AOA) and methanogenic archaea [23][24][25]. Thus, it is not very clear what the diversity and composition of the archaeal community as a whole is in the rhizosphere, especially under unique environmental stress such as that found on the Qinghai-Tibetan Plateau. In addition, the rhizosphere community structure is affected by the combination of environmental variables and interactions among microbial species [26,27]. However, given the unique cellular structure and speci c metabolic pathways of archaea that enable them to survive and even thrive under various adverse environments [12,28], several key questions about the archaeal community also need to be answered. One of the questions is about the assembly process of archaea in rhizosphere: is it governed by a deterministic process or stochastic process? The other question is how archaeal species interact with one another.
The Qinghai-Tibetan Plateau is the highest and largest plateau in the world, covering an area of approximately 2.5 × 10 6 square kilometres, and it has an average elevation of more than 4,000 metres [29]. Under the selection of extreme conditions in the Qinghai-Tibetan Plateau, most plants have evolved cold-tolerant characteristics [30], and soil microorganisms may form relatively distinctive microbial communities [31,32]. In this study, high-throughput sequencing of 16S rRNA gene amplicons was performed to exhaustively examine the archaeal communities derived from the rhizosphere of two native plants in the Qinghai-Tibetan Plateau. We aimed to address the following questions: 1) Are the structure of bulk soil and rhizosphere archaeal community different in the Qinghai-Tibetan Plateau? 2) What environmental factors signi cantly affect the archaeal community structure? 3) Whether the stochastic processes or the deterministic processes dominated the assembly of bulk soil and rhizosphere archaeal community? 4) Are the co-occurrence patterns of bulk soil and rhizosphere archaeal community different in the Qinghai-Tibetan Plateau?

Soil physicochemical properties
The soil physicochemical properties signi cantly differed between the rhizospheres of two plant species and the bulk soil (Table 1; Table S1). The pH varied from 7.84 to 7.91, and the lowest pH was in the bulk soil. The moisture of the two plant rhizospheres were similar and were lower than that of the bulk soil. The highest content of soil organic matter (SOM) was observed in the bulk soil, and a signi cant difference was detected only in the rhizosphere of P. crassifolia compared to the bulk soil (P < 0.05). In addition, there were no signi cant differences in the content of total nitrogen (TN), alkali-hydrolysable nitrogen (AN) and total phosphorus (TP) among the two plant rhizosphere and the bulk soil, but the content of ammonium nitrogen (NH 4 + -N) and available phosphorus (AP) in the two rhizosphere were signi cantly higher than they were in the bulk soil (P < 0.05).

Diversity And Community Composition Of Archaea
A total of 474,190 high-quality sequences were obtained with a median read count per sample of 39,516 (range: 30,420 − 54,538). The high-quality reads were clustered using > 97% sequence identity into 207 archaeal OTUs. The Good's coverage scores (in all cases above 99.9%) and the rarefaction curves showed clear asymptotes (Fig. S1), which together indicated a near-complete sampling of the archaeal community in this study.
The diversity indices of archaeal communities varied among the rhizospheres of two plant species and the bulk soil ( Table 2). The observed number of OTUs (Ob) was highest in bulk soil, followed by the rhizosphere of P. szechuanica, whereas the rhizosphere of P. crassifolia had lower numbers. Conversely, the Shannon index in the two plant rhizospheres were higher than they were in the bulk soil, and signi cant difference was identi ed only in the rhizosphere of P. szechuanica compared to the bulk soil (P < 0.05). The phylogenetic diversity (MNTD) of the two plant rhizospheres were similar, and their values were higher than that of the bulk soil. Principal coordinate analysis (PCoA) based on weighted UniFrac distances was performed to investigate the patterns of separation among archaeal microbiota. We clearly observed strong clustering of archaeal communities according to the different microhabitats (i.e., P. crassifolia, P. szechuanica rhizosphere and bulk soil). Moreover, the two plant rhizosphere samples were clearly distinguished from the bulk soil samples across the rst principal coordinate, while the separation between the rhizosphere of P. crassifolia and P. szechuanica was seen along the second principal coordinate, indicating that the largest source of variation in the archaeal communities is proximity to the root, followed by plant variety (Fig. 1A). Interestingly, PCoA analysis of βMNTD distances revealed that the largest source of variation is plant variety rather than proximity to the root (Fig. S2). Consistent with the result of PCoA analyses, ANOSIM analyses also revealed signi cant differences in the structure of archaeal communities among the rhizosphere of two plant species and the bulk soil (Table S2).
The relative abundance of archaeal OTUs at the phylum level was variable among the two plant rhizospheres and the bulk soil. The most dominant archaeal phyla across all samples were Thaumarchaeota, Unclassi ed_k_norank and Euryarchaeota, accounting for 92.46-98.01%, 1.35-6.01% and 0.56-1.18% of the pyrosequencing reads, respectively (Fig. 1B). Analysis of variance (ANOVA) showed signi cant enrichment of Thaumarchaeota in the rhizosphere microbiota of two plant species compared to that of the bulk soil (Dunnett test, P < 0.05). Conversely, the relative abundance of Unclassi ed_k_norank and Euryarchaeota in the rhizosphere microbiota of two plant species decreased but did not show signi cant differences compared with the abundance in the bulk soil (Table S3).
Moreover, LEfSe analysis was also performed to determine the taxa that most likely explains the variations among different samples. In the bulk soil, four groups of archaea were signi cantly enriched, namely, Thermoplasmata (the class, orders of Thermoplasmatales, and its family marine_Group_II to genus), unclassi ed_k_norank (from phylum to genus), norank_c_Soil_Crenarchaeotic_Group_SCG (from order to genus), group_C3 (from family to genus). In the P. crassifolia rhizosphere, a group of archaea was signi cantly enriched, namely, Thaumarchaeota (the phylum and its class soil_Crenarchaeotic). In the P. szechuanica rhizosphere, two groups of archaea were signi cantly enriched, namely, unclassi ed_c_Soil_Crenarchaeotic_Group_SCG (from order to genus), unknown_Order_c_Soil_Crenarchaeotic_Group_SCG (from order to genus) (Fig. 2).

Correlation Between Soil Properties And Archaeal Communities
Distance-based redundancy analysis (dbRDA) indicated the strong correlation between soil physicochemical characteristics and the structure of archaeal communities. The rst two axes of CAP could explain 27.12% and 13.43% of the total variation in archaea communities, respectively (Fig. 3). In line with the PCoA (weighted UniFrac) analysis, the rst axis (CAP1) could separate the rhizosphere samples from the bulk soil, and the second axis (CAP1) mainly distinguished the P. crassifolia rhizosphere from the P. szechuanica rhizosphere samples. The results of PERMANOVA analysis revealed that soil ammonium nitrogen (NH 4 + -N), soil organic matter (SOM) accounted for 35.1% and 28.5% of archaeal community differences, respectively, and niches (rhizosphere vs bulk soil) contributed 45.4% of the interpretation (Table S4). In addition, soil ammonium nitrogen (NH 4 + -N), available phosphorus (AP) and pH value were important environmental attributes signi cantly affecting the archaea community structure (Mantel test; r = 0.392, P = 0.026; r = 0.362, P = 0.030; r = 0.400, P = 0.028).
Further analyses revealed that soil properties had signi cant effects on the relative abundance of the archaea taxa at the class level. Soil pH value was positively correlated with the relative abundance of Unclassi ed_k_norank, Norank_p_Bathyarchaeota, and it was negatively correlated with Soil_Crenarchaeotic_Group_SCG, Methanobacteria. Ammonium-nitrogen (NH 4 + -N) was positively correlated with the relative abundance of Thermoplasmata. Soil total phosphorus (TP) was positively correlated with the relative abundance of Methanobacteria. Available phosphorus (AP) was negatively correlated with Unclassi ed_k_norank (Fig. 4).
Assembly processes of archaeal microbiota in rhizosphere and bulk soil The phylogenetic tree of archaea recovered from all samples were relatively well classi ed according to the major lineages, and the local support values on the branches were relatively high ( Fig. S3), suggesting the archaeal phylogenetic tree was reliable. Additionally, the phylogenetic signal showed that there was a signi cant relationship between ecological similarity and phylogenetic relatedness (Fig. S4). Thus, it was feasible to gain insight into the assembly processes of archaea microbiota based on phylogenetic analysis [33].
We clearly observed that the NTI values of archaea microbiota from all samples were less than − 2, in which the lowest mean NTI value was detected in the bulk soil (Fig. 5A), suggesting that archaeal communities were phylogenetically over-dispersed, especially in the bulk soil, and it also suggested that deterministic processes mainly regulate the assembly of archaeal communities. In addition, the lowest mean βNTI value for archaea was found in the bulk soil, and it was signi cantly lower than zero (66.67% of βNTI values ranging between 0 and − 2, 33.33% of βNTI values less than − 2), suggesting the phylogenetic turnover was less than what would be expected by chance. The mean βNTI value in the P. crassifolia rhizosphere was signi cantly higher than zero, suggesting the phylogenetic turnover was higher than what would be expected by chance. However, the mean βNTI value in the P. szechuanica rhizosphere was not signi cantly different from zero (Fig. 5B). These results further indicated that deterministic processes play a stronger role in the phylogenetic turnover than stochastic processes. Additionally, the species rank abundance distribution models also revealed that archaea communities in all samples were followed 'niche theory' models (Table 3).

Co-occurrence Network Structure Of Archaea Microbiota
Three co-occurrence networks were constructed for all sample types (Bulk soil, P. crassifolia and P. szechuanica rhizosphere) to illustrate potential biotic interactions among archaea taxa (Fig. 6). All the networks were signi cantly different from the random networks with the identical numbers of nodes and edges (Table S5), suggesting that the network structures were non-random and reliable.
We found that the networks in the rhizosphere of two plant species were obviously different from that in the bulk soil ( Fig. 6; Table S5). The number of edges, average degree and average clustering coe cient of the networks in two rhizosphere were lower than they were in bulk soil, indicating that rhizosphere assemblages of two plant species formed lower complex archaea networks compared with that of the bulk soil. The ratio of negatively correlated edges between OTUs in the P. crassifolia rhizosphere (30.8%) and the P. szechuanica rhizosphere (27.6%) were profoundly higher than that of in bulk soil (20.1%), which could be interpreted as increased competitions among archaea taxa in the rhizosphere environment. We also observed a high proportion of unclassi ed_k_norank in all networks (Bulk soil, P. crassifolia and P. szechuanica rhizosphere), accounting for 52.0%, 52.0% and 62.0%, respectively. In addition, the majority of unclassi ed_k_norank were highly connected in the networks. Thus, it could be inferred that unclassi ed_k_norank is very crucial for the stability of archaea network structures in all samples.

Discussion
Variation of archaea community structures between rhizosphere and bulk soil In our study, the archaeal community structures were signi cantly different between the rhizospheres of two plant species and the bulk soil according to both taxonomic and phylogenetic metrics ( Table 2; Table  S2; Fig. 1; Fig. S1). These differences could be explained in two ways. First, plant roots could release a variety of carbon exudates including sugars, amino acids, organic acids, mucilage and root border cells [34,35]. These exudates are available nutrients and energy for microbial activities [36], making the microbial community structures in the rhizosphere differed from what is found in the bulk soil [5]. In particular, the differences of archaeal communities observed between P. crassifolia and P. szechuanica rhizosphere could be attributed to the differences in carbon exudates released by two plant species [37].
This effect of plant root exudates is also re ected in bacterial and fungi community in rhizosphere [38][39][40][41]. A second possible explanation is that the rhizosphere of plants could form oxygen-depleted microniches owing to root respiration, which provide a special habitat for some important archaea taxa such as methanogens and ammonium-oxidizing archaea [25,42].
Analysis of archaeal community composition revealed that the archaeal communities were dominated by Thaumarchaeota, accounting for 92.46-98.01% of sequences in this study ( Fig. 1B; Table S3). This nding was in agreement with previous ndings from the research also conducted in the Qinghai-Tibetan Plateau, which showed that the dominant archaeal phylum was Thaumarchaeota, accounting for 79.27% of sequences [31] . Thaumarchaeota have been detected in a variety of habitats [43][44][45], and identi ed as a novel archaeal phylum in 2008 [43]. Many studies have suggested that Thaumarchaeota species possess ammonia oxidizing abilities and are considered to play an important role in nitrogen cycling [46,47]. In our study, the relative abundance of Thaumarchaeota in the two different plant rhizosphere were signi cantly higher than they were in the bulk soil ( Fig. 1B; Table S3), and all the eight biomarkers in rhizosphere by the LEfSe analysis also belong to the phylum Thaumarchaeota (Fig. 2). These ndings collectively indicated that the nitrogen metabolism activities occurred in the rhizosphere might be higher than that in that in the bulk soil.

Important Drivers Of Archaea Communities
Combined with the analysis of dbRDA and Mantel test showed that soil ammonium nitrogen (NH 4 + -N), available phosphorus (AP) and pH were signi cantly correlated with the archaeal community structures (Fig. 3). This observation in NH 4 + -N agreed with previously reported results in the study by Norman and Barrett [48], which documented NH 4 + -N as a metabolic substrate that drives the distribution patterns in richness of ammonia-oxidizing archaea (AOA). Moreover, our results revealed that the content of NH 4 + -N in the rhizosphere was signi cantly higher than it was in the bulk soil (Table 1), and the increase of NH 4 + -N concentration in the rhizosphere may be due to the enrichment of diazotrophic bacteria in the rhizosphere, which can convert atmospheric N 2 into ammonium via biological nitrogen xation [49,50], although the results of this study cannot be directly con rmed. The relative abundance of the class Thermoplasmata was positively correlated with soil NH 4 + -N (Fig. 4). These results corroborate the opinion that NH 4 + -N may play a signi cant role in shaping the archaeal community structures. Numerous studies have suggested that soil pH is a major driver of the community structure of bacteria, fungi or diazotrophs [51][52][53], but there seems to be no consensus on archaea [48,54]. This contradiction could be explained by a plausible interpretation: soil pH indirectly affects the abundance of major archaeal taxa mainly by regulating the availability of substrates such as NH 4 + , CO 2 , and CH 3 COOH [48,55,56], so the correlation may vary with samples or by region. The signi cant correlation in our study could be attributed to more available NH 4 + -N regulated by pH to the phylum Thermoplasmata [57]. Soil available phosphorus (AP) has been reported to be a limiting factor for the growth of plants or microorganisms [7,58]; thus, it may directly or indirectly affect archaeal communities.

Deterministic Processes Govern The Assembly Of Archaeal Communities
The phylogenetic metrics (NTI and βNTI) signi cantly differed from a null distribution in all samples (Fig. 5), indicating the dominant role of deterministic processes, according to the framework in the study by Stegen et al. [33]. Moreover, this interpretation was supported by the results of the species rank abundance distribution models in our study ( Table 3). The mean NTI values were signi cantly lower than zero in all samples (Fig. 5A) and provided concrete evidence that the archaeal communities were more phylogenetically over-dispersed than expected as a result of chance [59,60]. Previous studies have shown that the competition among species would become more frequent where there was greater niche similarity and would subsequently lead to the coexistence of distantly phylogenetically related species [61,62]. In our study, most of the soil variables (except for NH 4 + -N and AP) were similar (Table 1), and microorganisms competing strongly for nutrients or water in the Qinghai-Tibet Plateau suffered from its low temperature and strong ultraviolet radiation [63,64], which was supported by the results of cooccurrences ( Fig. 6; Table S5). These factors may explain why the archaeal communities were phylogenetically over-dispersed. Furthermore, the NTI values in the rhizosphere were greater than that in the bulk soil but not signi cant, which might indicate that stochastic processes may still play a minor role [65]. We also noticed that the rhizosphere βNTI values were signi cantly greater than those measured in the bulk soil (Fig. 5B), suggesting that the phylogenetic turnover of archaea in the rhizosphere were higher than what was in the bulk soil [50,66]. This could be attributed to dynamic rhizosphere microhabitats potentially stimulating the activities and evolutions of archaeal species [3].

Distinct Archaeal Networks In Rhizospheres And Bulk Soil
Previous studies have found that bacterial or fungal networks in the rhizosphere were more [40,67] or less [27,68] complex than what was found in bulk soil. In the present study, we found that the cooccurrence networks of archaea in the rhizosphere of two plant species had fewer correlations when compared with that of the bulk soil ( Fig. 6; Table S5), which could also be explained through two plausible interpretations. First, the archaea possess distinctive metabolic pathways and enzymes that enable them to survive and thrive under extreme or nutrient-poor environments [12,13]. Thus, the archaeal species might have a lower nutritional dependence on root exudates than bacteria or fungi. On the other hand, the rhizosphere bacterial and fungal species are likely to accelerate the consumption of substrates required by archaea, and even the plants themselves may be competitors for microorganisms under environmental stress [69,70], thus reducing the interactions or niche sharing among archaea. This interpretation was also supported by the nding that higher numbers of negative links occurred in the rhizosphere networks ( Fig. 6; Table S5). Collectively, these two mechanisms may introduce fewer trophic interactions among archaea in the rhizosphere than that in the bulk soil. Our results also revealed that as an unclassi ed archaeal group, unclassi ed_k_norank occupies a high proportion (Fig. 6) in all networks despite its low abundance in the community composition ( Fig. 1B; Table S3). Even though we do not yet know the speci c ecological functions of this unclassi ed archaeal group, it is clear that it may play an important role in maintaining the stability of community structure and function [71].

Conclusions
In summary, this study provides insight into the structure, assembly and co-occurrence patterns of the rhizosphere archaeal communities in the Qinghai-Tibetan Plateau. The results showed that archaeal community structures in the rhizosphere of two plant species signi cantly differed from that in the bulk soil. Soil ammonium-nitrogen (NH 4 + -N), soil organic matter (SOM), available phosphorus (AP) and pH were important drivers of the archaeal communities. Deterministic processes dominated the assembly of archaeal communities across all samples. The network structures of the archaeal community in the rhizosphere were less complex than they were in the bulk soil. We also identi ed an unclassi ed archaeal group (unclassi ed_k_norank) that may be crucial for the interrelationships among archaeal species. Future research should further investigate the interaction between archaea and other microorganisms such as bacteria, fungi and protists in the rhizosphere, and work to understand the role of archaea in plant survival and growth under low-temperature stress.

Site and sampling
A trees eld trial located in the northeast portion of the Qinghai-Tibetan Plateau (31°32'N, 92°00'E, 4531 m a. s. l), which has a plateau sub-frigid monsoon semi-arid climate with an average annual temperature of -2.2 °C and a mean annual precipitation of 458 mm. This eld trial was established in April 2010 and contains two native alpine tree species (Picea crassifolia and Populus szechuanica var. tibetica).
In July 2017, we selected the surviving and well-growing trees (P. crassifolia about 2.5 m, P. szechuanica about 4.5 m) for sample collection, and repeated four trees for each plant species. In order to ensure the representativeness of the samples, three subsamples of ne roots (< 2 mm) were carefully collected from different positions in the rhizosphere of each selected tree at the depth of 5-15 cm below ground level. The homogeneous rhizosphere soil was obtained from the combined ne root samples of each tree according to the procedure described in a previous study [38]. The bulk soil was collected from four treeless quadrats (3 m × 3 m) at the depth of 5-15 cm below ground level. Each quadrat is about 10 m away from the sampled trees, in which ve soil subsamples were obtained and combined into a representative bulk soil sample. All soil samples were handpicked to remove roots and impurities, and then divided into two subsamples. One portion was air dried and sieved through 2 mm meshes for soil property analyses, and the other portion was stored at -80 °C for DNA extraction.

Soil Physicochemical Properties Analysis
Soil physicochemical properties in both rhizosphere and bulk soils were analysed according to Zhou et al. [72]. Brie y, soil moisture was quanti ed gravimetrically by drying fresh soils in 105 °C for 48 h. Soil pH was measured by a pH meter in a soil suspension with an air-dried soil to water radio of 1: 2.5 mass/volume. Soil organic matter (SOM) was determined by the potassium dichromate oxidation titration method. Total nitrogen (TN) was measured by the Kjeldahl digestion method. Total phosphorus (TP) was determined by the Mo-Sb anti-spectrophotometric method. Alkali-hydrolysable nitrogen (AN) was measured by the alkali-hydrolysed diffusing method. Ammonium nitrogen (NH 4 + -N) was measured using indophenol blue spectrophotometry. Available phosphorus (AP) was extracted with a NH 4 F/HCl solution, which was then determined using a UV-visible spectrophotometer.
Dna Extraction, Pcr Ampli cation, And Sequencing

Sequence Processing
Raw sequences yielded from Illumina sequencing were processed using QIIME 1.9.1 [74]. Paired-end reads were joined with fastq-join, demultiplexed and quality ltered with default parameters [75]. Brie y, sequences with a quality score < 20 or with any truncated reads shorter than 50 bp were removed.
Operational taxonomic units (OTUs) were clustered with 97% similarity cutoff using UPARSE 7.1 and chimeric sequences were identi ed and removed using UCHIME. The taxonomy of each 16S rRNA gene sequence was analysed by an RDP Classi er algorithm (http://rdp.cme.msu.edu/) against the Silva (SSU123) 16S rRNA database using a con dence threshold of 80%.

Data analysis
All statistical analyses were carried out using R (v 3.5.1, The R Core Team, 2018) unless stated otherwise. α-diversity in each sample was calculated as the observed number of OTUs (Ob), the Shannon diversity and the phylogenetic diversity (MNTD) indices. Signi cant differences in the variance of α-diversity and microbial abundance data were examined using one-way analysis of variance, and post hoc comparisons were conducted by the Dunnett test at the 5% level. The differences in archaeal community composition based on Weighted UniFrac and βMNTD distances were illustrated with PCoA ordination plots using the 'cmdscale' function from the vegan package. To statistically support the archaeal clustering patterns resulted from PCoA analysis, different samples were compared by ANOSIM analysis using the vegan package. Additionally, we performed linear discriminant analysis (LDA) coupled with effect size measurements (LEfSe) analysis to investigate statistically representative biomarkers between different samples.
Distance-based redundancy analysis (dbRDA) was carried out using the 'rda' function from the vegan package to explore the relationships between soil physicochemical properties and archaeal community composition. Furthermore, associations between soil properties and nine archaeal classes were evaluated by Pearson correlation analysis at the 5% level.
To evaluate the assembly processes of the archaeal community, the phylogenetic signal of each sample was rst tested by following the procedure described by a previous study [65]. Brie y, environmental optima for each OTU with respect to all physicochemical variables were calculated. The correlation coe cients between phylogenetic distances and differences in environmental optima were calculated by phylogenetic Mantel correlograms [50], and the signi cance of these correlations were examined with 999 randomizations using the 'mantel.correlog' function from the vegan package. The phylogenetic diversity within each sample was calculated as the mean nearest taxon distance (MNTD) and nearest taxon index (NTI) using the 'mntd' and 'ses.mntd' functions from the picante package [76] Note that MNTD refers to the phylogenetic distance between each OTU and its closest relative also found per sample, and NTI measures the deviation of observed MNTD from MNTD in a null model with 999 randomizations. For NTI > + 2 (NTI < − 2) in a single community or a mean NTI > 0 (NTI < 0) signi cantly across all communities indicates coexisting taxa are more closely related (phylogenetic clustering) or more distantly related (phylogenetic over-dispersion) than can be expected by chance [33]. The pairwise phylogenetic turnover between communities was calculated as βMNTD and βNTI using the 'comdistnt' function from the vegan package [76]. βNTI > + 2 (βNTI < − 2) between one pair of communities or mean βNTI > 0 (βNTI < 0) signi cantly in all pairs of communities indicates greater (or less) than expected phylogenetic turnover, respectively [65]. If the observed βMNTD values does not signi cantly deviate from the null βMNTD distribution [66], it suggests that stochastic processes predominate phylogenetic community composition. In addition, to verify the results from phylogenetic analyses, ve models representing niche theory (Break Stick, Pre-emption, Lognormal, Zipf, Zipf-Mandelbrot) and ZSM representing neutral theory were performed using the function 'rad t' from the R package 'vegan' or TeTame [60]. Akaike Information Criterion (AIC) was used to evaluate the tting quality of each statistical model, where the lower AIC value indicated a better t for the model [77].
Network analyses based on Spearman's rank analysis were carried out with the 'WGCNA' package [78], and structural attributes of the overall networks including average degree, clustering coe cient and average path distance were calculated in the 'igraph' package. The 50 most abundant OTUs of the archaea community in each sample were selected, and the co-occurrence patterns of archaea communities were explored based on strong and signi cant correlation (Spearman's |r| > 0.6, P < 0.05). Finally, the constructed networks were visualized using Gephi 0.9.2 [79].

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
All data generated or analyzed during this study are included in this published article and its supplementary information les.

Competing interests
The authors have no con icts of interest with any parties or individuals.    Ordination plots of the results from distance-based redundancy analysis (dbRDA) to explore the correlation between soil properties and archaeal community structure.

Figure 4
Heat map showing the Pearson correlation between soil properties and the relative abundance of the archaea taxa at the class level. The left side of the legend is the colour range of different R-values. The value of P < 0.05 and P < 0.01 is marked with "*" and "**", respectively. The co-occurrence network of archaea communities based on correlation analysis. Each node represents an individual OTU coloured by taxonomy at phylum level, and the size of each node is proportional to the degree. The connection stands for a strong and signi cant (Spearman's |r| > 0.6, P < 0.05) correlation. Red edges represent positive, blue edges represent negative correlation.