Soil physicochemical properties
The soil physicochemical properties significantly 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 significant difference was detected only in the rhizosphere of P. crassifolia compared to the bulk soil (P < 0.05). In addition, there were no significant 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 (NH4+-N) and available phosphorus (AP) in the two rhizosphere were significantly 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 significant difference was identified 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 first 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 significant 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, Unclassified_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 significant 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 Unclassified_k_norank and Euryarchaeota in the rhizosphere microbiota of two plant species decreased but did not show significant 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 significantly enriched, namely, Thermoplasmata (the class, orders of Thermoplasmatales, and its family marine_Group_II to genus), unclassified_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 significantly enriched, namely, Thaumarchaeota (the phylum and its class soil_Crenarchaeotic). In the P. szechuanica rhizosphere, two groups of archaea were significantly enriched, namely, unclassified_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 first 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 first 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 (NH4+-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 (NH4+-N), available phosphorus (AP) and pH value were important environmental attributes significantly 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 significant effects on the relative abundance of the archaea taxa at the class level. Soil pH value was positively correlated with the relative abundance of Unclassified_k_norank, Norank_p_Bathyarchaeota, and it was negatively correlated with Soil_Crenarchaeotic_Group_SCG, Methanobacteria. Ammonium-nitrogen (NH4+-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 Unclassified_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 classified 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 significant relationship between ecological similarity and phylogenetic relatedness across short phylogenetic distances (Fig. S4). Thus, we calculated NTI and βNTI because both of these metrics emphasize phylogenetic relationships across short phylogenetic distances .
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 significantly 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 significantly 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 significantly 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 significantly 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 tree species were obviously different from that in the bulk soil (Fig. 6; Table S5). The number of edges, average degree and average clustering coefficient 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 unclassified_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 unclassified_k_norank were highly connected in the networks. Thus, it could be inferred that unclassified_k_norank is very crucial for the stability of archaea network structures in all samples.