Soil sampling and the determination of soil physicochemical properties were conducted as previously described (Liu et al., 2014). In brief, based on the China Black Soil Ecology database (http://www.blackland.csdb.cn), 26 arable soil samples were collected across the black soil zone of northeast China in September 2012 (Figure S1). Each soil sample was a composite of 10 soil cores that were randomly collected from tillage layers (0–20 cm) within an area of ~100 m². DNA was extracted from 0.5 g moist soil using an E.Z.N.A Soil DNA kit (OMEGA, USA) according to the manufacturer's instructions. The extracted DNA was diluted in TE buffer (10 mM Tris-HCl, 1 mM EDTA, adjust pH to 7.0) and stored at−20°C until used for additional analyses. Soil physicochemical properties are listed in Table S1.

Archaeal abundance was quantified using real-time PCR targeting the hypervariable V3–V4 regions of the 16S rRNA gene with the primers Arch519f (5′-CAG CCG CCG CGG TAA-3′) (Øvreås et al., 1997) and Arch915r (5′-GTG CTC CCC CGC CAA TTC CT-3′) (Stahl and Amann, 1991; Inceoglu et al., 2015). Each PCR reaction contained 10 μL of SYBR Premix Ex Taq^TM (Takara, Dalian, China), 0.4 μL of 10 μM each forward and reverse primers, 2 μL of extracted soil DNA, and the volume was reached to 20 μL with sterilized MilliQ water. Amplification was performed in a LightCycler^®480 (Roche Applied Science) using initial denaturation at 95°C for 30 s, followed by 30 amplification cycles (95°C for 30 s, 60°C for 40 s, 70°C for 30 s), and a final extension at 50°C for 5 min for cooling. Standard curves were generated using a 10-fold dilution series of known concentrations using a plasmid containing the archaeal 16S rRNA gene amplicon. The copy number of the archaeal genes was calculated by a regression equation to convert the cycle threshold value to the known copy numbers in the standard curves. The copy numbers of the archaeal genes were converted into per gram of dry soil as described in our previous study (Liu et al., 2018). In addition, the ratio of between the archaeal and bacterial 16S rRNA gene copies was also calculated using abundance data of bacterial 16S rRNA genes reported previously (Liu et al., 2016a).

Archaeal 16S rRNA genes were PCR amplified using the primer pairs Arch519f/Arch915r (Bai et al., 2017), with the forward primer being modified to contain a unique 8 nt barcode at the 5′ end. PCR reactions were performed in triplicate in a 25 μL mixture containing 1.0 μL purified DNA template (1–10 ng), 0.5 μL of each primer (10 μM), and 23 μL of Platinum PCR SuperMix (TransGen Biotech Co. Ltd., Beijing, China). Each sample was amplified under the following conditions: initial denaturation at 95°C for 5 min, followed by 35 cycles (95°C for 1 min, 60°C for 45 s, 70°C for 1 min) with a final extension at 72°C for 5 min. Each PCR product was pooled and purified using an Agarose Gel DNA purification kit (TaKaRa, Dalian, China). Equimolar amount of the PCR products were combined into one pooled sample and prepared for sequencing (2 × 300) using an Illumina MiSeq platform at Majorbio BioPharm Technology Co., Ltd. (Shanghai, China).

The raw data were processed and analyzed using QIIME Pipeline (http://qiime.org/tutorials/tutorial.html) (Caporaso et al., 2010b). Original paired reads were joined with FLASH (fast length adjustment of short reads) software (Magoč and Salzberg, 2011). Low quality sequences which shorter than 200 bp, and with an average base quality score of <20 were removed for further analysis. The trimmed sequences were chimera-detected and removed with the flags - -non_chimeras_rentention = intersection basing on a combination of de novo and reference-based chimera checking. The high-quality sequences were clustered into OTUs with UCLUST based on a 97% similarity level using QIIME's pick_ open_reference_otus.py (Edgar et al., 2011). The taxonomic classification of each phylotype (OTU) was aligned using the Python Nearest Alignment Space Termination (PyNAST) tool (Caporaso et al., 2010a) and a relaxed neighbor-joining tree built by FastTree (Price et al., 2009). The taxonomic identity of each phylotype was determined using the Ribosomal Database Project (RDP) classifier with a confidence threshold of 0.80 (http://rdp.cme.msu.edu/) (Cole et al., 2005). After the taxonomies had been assigned, the OTUs not matched to the archaea were removed from the dataset. The raw reads obtained in this study were submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) with accession number SRP139400.

The OTU richness was used to compare soil archaeal alpha diversity. The SPSS software (version 18.0 for Windows) was adopted to conduct Pearson's correlation analysis to identify possible correlations between abundance of archaeal genes, alpha diversity, relative abundance of the taxonomic subgroups (relative abundance >0.3% at least in one sample) and soil characteristics. In order to clarify the taxonomic status of Thaumarchaeota and archaea, the phylogenetic identities of each phylotype obtained in this study were aligned with the taxonomy-determined sequences obtained from the NCBI database using ClustalX 1.81 (Thompson et al., 1997). Briefly, the alignment analysis was performed with ClustalX 1.81, and a rooted neighbor-joining tree of Thaumarchaeota and an unrooted phylogenetic tree of archaea was constructed using Molecular Evolutionary Genetic Analysis software (MEGA 7.0) with 1,000-fold bootstrap support (Kumar et al., 2016). Meanwhile, an unrooted phylogenetic tree was illustrated by the interactive Tree of Life online program (Letunic and Bork, 2006).

To correct for sampling effort, we used a randomly selected subset of 8,300 archaeal sequences per sample to calculate the archaeal community compositions. Non-metric multidimensional scaling (NMDS) was used to analyze beta diversity in different soil samples, which was conducted on the basis of the calculated weighted UniFrac distances (Lozupone and Knight, 2005). Canonical correspondence analysis (CCA) was used to test the significant effects of environmental variables that drive the difference in the archaeal communities with latitude, soil water content, total C, N, and P, C/N, pH, available K and P, NH4+-N, and NO3--N as the explanatory factors. Only soil factors which had significant effects on the archaeal communities were retained for the CCA analysis that was selected using the bioenv () function in the R environment. In addition, a VPA was conducted to quantify the relative contributions of the special factors and environmental variables on the archaeal community compositions using the “vegan” package in R 2.3.5 (R Development Core Team, 2010).

The archaeal abundance of 16S rRNA gene in all soil samples ranged from 4.04 × 10⁷ to 26.18 × 10⁷ copies per gram of dry soil (Table S1). The ratio between the archaeal and bacterial 16S rRNA gene copies ranged from 0.032 to 0.085. Pairwise correlation analysis shown that the abundance of archaea highly correlated with soil pH (r = 0.516, p = 0.007), soil C content (r = 0.821, p < 0.0001), soil N content (r = 0.757, p < 0.0001), soil water content (r = 0.676, p < 0.0001), available K content (r = 0.609, p = 0.001), and soil P content (r = 0.629, p = 0.001) (Figure S2).

In total, 735,020 high-quality archaeal sequences were obtained from all 26 samples. Of these, 74.12% could be classified as archaeal sequences by classifier alignment using the RDP database, and 8,313–38,342 sequences were obtained per sample (mean 25,885). The read lengths ranged from 381 to 419 bp with a mean of 384 bp.

Thaumarchaeota was the dominant archaeal phylum accounting for an average of 94.15% of the total observed sequences (ranging from 85.44 to 98.86%). Euryarchaeota was the second most abundant phylum with an average of 5.72% across all the samples (ranging from 1.14 to 13.29%). Crenarchaeota was sporadically detected in 10 soil samples with a relative abundance below 0.31%, with the exception of a location of LS with a relative abundance of 2.74% (Table S2). None of the other archaeal phyla were detected in the black soils.

Within the phylum of Crenarchaeota, only two orders Desulfurococcales and Thermoproteales were occasionally observed in the black soils (Table S2). Within the Euryarchaeota, 10 orders were observed with Methanomassiliicoccales detected in all soil samples with a relative abundance ranging from 0.75 to 5.87% (average of 3.03%), followed by Methanobacteriales ranging from 0.24 to 5.50% (average of 2.04%). Other orders were sporadically detected in some soil samples with very low abundances. Within Thaumarchaeota, three orders were detected, of which Nitrososphaerales was the most abundant with a relative abundance that ranged from 58.52 to 97.55% (average abundance of 87.35%), followed by Thaumarchaeota_norank with an average of 6.27%; Nitrosopumilales was only detected in 7 of the 26 samples with a relative abundance below 0.6%, with the exception of two samples of CT1 and LS with a relative abundance of 5.41 and 7.81%, respectively (Table S2).

In total, 185 archaeal OTUs were detected in this study, and the OTU number ranged from 41 to 130 (Table S1). Among the 185 OTUs, 3, 31, and 151 were found within the phyla Crenarchaeota, Thaumarchaeota and Euryarchaeota, respectively. Figure 1 shows the phylogenetic tree of the OTUs observed in this study, which indicated that the diversity of Euryarchaeota was the highest archaeal phylum in the black soils. The relative abundances, taxonomy information and represented sequences of all the OTUs are shown in Table S3. The closest relative of each OTU identified by BLAST search are shown in Table S4.

The phylogenetic positions of the 31 OTUs with the sequence numbers within Thaumarchaeota are shown in Figure S3. Two, 2, 8, and 19 OTUs fell into group 1.1a, group 1.1a-associated, group 1.1c, and group 1.1b, respectively. The sequences belonging to group 1.1b were the most abundant, accounting for 92.35% of the sequences of Thaumarchaeota, followed by group 1.1a-associated (6.93%), group 1.1a (0.59%), and group 1.1c (0.13%). Seven OTUs (OTU1~OTU6, OTU1120) were abundant members with total sequences accounting for 80.33% of the Thaumarchaeota sequences.

Pairwise correlation analysis showed that the alpha-diversity (OTU number) of the archaea had no relationship with any soil parameters. However, the OTU number within the phyla Thaumarchaeota and Euryarchaeota showed negative (r = −0.442, p = 0.024) and positive (r = 0.406, p = 0.040) correlations with the soil pH, respectively (Figure 2).

The relationship between the sampling latitude, soil properties and relative abundance of archaea at the different taxonomy levels was examined using pairwise correlation analysis (Table 1). None of the three phyla had a significant relationship with any environmental variables detected. At the order, genus and OTU levels, only the members belonging to Thaumarchaeota had significant correlations with some environmental variables. For example, the relative abundance of Nitrososphaerales was positively correlated with the sampling latitude, soil pH, TC, and TN concentrations, while the relative abundance of Thaumarchaeota_norank exhibited the reverse tendency. At the genus level, group 1.1a-associated and group 1.1 showed an opposite tendency in their change with the sampling latitude, soil pH, total C and total N concentrations. At the OTU level, 17 OTUs were significantly correlated with the environmental variables. Among these, the abundant members of OTU1 and OTU5 correlated positively, while OTU3 correlated negatively with the sampling latitude; OTU1120, OTU2, and OTU3 were negatively correlated, while OTU6 positively correlated with the soil pH; in addition, OTU1 also showed a significant positive relationship with the total C, total N and soil moisture concentrations, while OTU3 showed the opposite trend (Table 1).

The archaeal community structures were illustrated using weighed NMDS plots based on the gradient of the soil pH and the soil total C content (Figure 3). The NMDS1 and NMDS2 axis score was negatively correlated with the soil pH (r = −0.836, p < 0.001) and positively correlated with the soil total C (r = 0.649, p < 0.001), respectively (Table S5). The CCA analysis indicated that the archaeal community structures differed among the sampling locations (Figure 4). CCA1 and CCA2 explained 28.01 and 16.81% of the community structures, respectively. Among the environmental factors examined, the soil pH was the dominant factor in shifting community structures along the CCA1 axis, while the latitude of the sampling locations was the dominant factor driving the changes of communities along the CCA2 axis. In addition, a VPA analysis showed that the geographic distances explained 19.01% of the community variation, and the selected soil parameters explained 32.35% with soil pH being the master variable explaining 20.89% of the community variation, leaving 47.55% of the variation unexplained.

In this study, compared with the data of bacterial abundance reported previously (Liu et al., 2016a), the ratio of archaeal to bacterial 16S rRNA gene copies ranged from 0.032 to 0.085, indicating that the archaeal abundance accounted for 3.11–7.80% of the total prokaryotes in the black soils. The abundance of the archaea in the agricultural soils varied between studies. Bates et al. (2011) used universal primers 515f/806r for HTS sequencing and found that the archaea comprised 2% of all the sequences with a range from 0 to >10% in individual soils. Bengtson et al. (2012) reported that the ratio of archaeal to bacterial 16S rRNA gene copy numbers ranged from 0.02% at low pH values to more than 7% at pH 8 in the Hoosfield acid strip at Rothamsted Research, UK. Cao et al. (2012) also reported that archaea accounted for 0.20–9.26% of the total prokaryotes in the selected Chinese soils. Therefore, the archaeal abundance across black soil zoon of northeast China is consistent with above-mentioned studies.

The influence of soil pH on the abundance of the archaea varied with the studies. Tripathi et al. (2013) detected the archaeal abundance to be negatively correlated with the soil pH across 27 soil samples collected from tropical forest and non-forest sites. However, Bengtson et al. (2012) found a strong positive relationship between the archaeal abundance and the soil pH in the temperate region. The soil samples used for this study were collected from the temperate region, and our finding is consistent with that of Bengtson et al. (2012). In addition to the soil pH, the soil organic C and total N were positively correlated with the archaeal abundance in several selected Chinese soils (Cao et al., 2012), which partially agreed with the findings of this study, showing that the soil C, N, and P concentrations, available K concentration and soil moisture content were highly positively correlated with the archaeal abundance (Figure S2). This finding suggests that there are multiple soil factors controlling the population of archaea in the black soils.

Thaumarchaeota was the dominating archaeal phylum in the black soils, which is consistent with the findings of Tripathi et al. (2015) who investigated 24 samples of topsoil collected from a tropical biome and 25 soils from a temperate biome. Our results are also consistent with the findings of Shi et al. (2016) who found that Thaumarchaeota was the most abundant phylum in the soils of eastern Tibetan Plateau. However, at the lower taxonomic level, Thaumarchaeota group 1.1b was the most abundant in the black soils (average of 87.23% of the total archaeal sequences), while in the results of Tripathi et al. (2015), Thaumarchaeota group 1.1b accounted for 48.7% of all the archaeal sequences, but the proportion of group 1.1b in Tibetan Plateau was not reported (Shi et al., 2016). This finding indicated that the members of Thaumarchaeota group 1.1b were the overwhelming archaea in the arable black soils. Recently, Lu et al. (2017) showed that the relative abundance of Thaumarchaeota group 1.1b was higher in the agricultural soils than in the forest soils.

Considering group 1.1b of Thaumarchaeota was taxonomically classified as Nitrososphaera (Bomberg, 2016), the relative abundance of Nitrososphaera (group 1.1b) in this study was comparable to the data observed in our former study of the AOA, in which the Nitrososphaera cluster and the Nitrososphaera sister cluster together comprised 89.34% of the AOA sequences across the same soil samples (Liu et al., 2018). Thus, we hypothesized that the majority of archaea in black soils contribute to nitrification.

Thaumarchaeota group 1.1a and group 1.1a-associated were phylogenetically classified into Nitrosopumilus and Nitrosotalea, respectively (Table S2; Bomberg, 2016). The findings of this study and our previous study of the AOA (Liu et al., 2018) indicated that those groups were minor, and the relative abundances of those groups are comparable between the two studies. For example, the relative abundance of the Nitrosospumilus cluster and Nitrosotalea comprised averages of 1.98 and 8.68% of the AOA sequences, respectively (Liu et al., 2018). Similarly, Thaumarchaeota groups 1.1a and 1.1a-associated accounted for averages of 0.54 and 6.23% of the archaeal sequences in this study (Table S2). The consistent results between the two studies suggest that the results observed in this study were correct and that the primer set used was suitable for analysis of the archaeal communities at least in the black soils.

Thaumarchaeota group 1.1c was frequently detected in acidic agricultural soils (Lehtovirta et al., 2009) and acidic forest soils (Jurgens et al., 1997; Bomberg and Timonen, 2009; Weber et al., 2015). However, in this study, this group only accounted for an average of 0.12% of the total archaeal sequences (Table S2), suggesting that group 1.1c is a rare member of the archaea in the arable black soils, even though the pH of most of the soils used for this study was below 6.0 (Table S1). Unlike Thaumarchaeota groups 1.1a and 1.1b, the phylogenetic position of group 1.1c was unclear, since no representative of this group has yet been obtained in pure culture (Weber et al., 2015; Bomberg, 2016). A recent study showed that group 1.1c are heterotrophic archaea, and its growth was not related to ammonia oxidation (Weber et al., 2015).

Although the relative abundance of Euryarchaeota only comprised 5.75% of the total sequences (Table S2), the OTU number within this phylum was greater than that in the Thaumarchaeota (Figure 1), suggesting that the diversity of Euryarchaeota is at the highest level among the archaeal phylum in the black soils. An interesting finding within the Euryarchaeota is that the majority of sequences are related to the methanogenic archaea (Table S2). Methanogenic archaea are usually considered to exist in highly reduced and anoxic conditions, such as wetlands, paddy fields, marine, or lake sediments, and the intestinal tracts of human and animals (Garcia et al., 2000; Liu and Whitman, 2008). However, the samples used in this study were collected from upland arable soils, where methane emission is not thought to be an important ecological process. We hypothesize that the diversity of the methanogenic archaea is related to the anaerobic conditions in localized micro sites, anaerobic seasonal conditions, or due to land conversion. A large number of OTUs grouped into the methanogenic Euryarchaeota were also found in upland soils across from South and North China (Hu et al., 2013).

In contrast to the Euryarchaeota, only 31 OTUs of Thaumarchaeota were observed across the black soils. Among these, seven OTUs were overwhelmingly dominant (Table S3), which suggested that the dominating members of archaea in the black soils are simplified. This finding is consistent with the observations of Bates et al. (2011) that two OTUs (DSC1 and DSC2) were highly abundant across 146 soils obtained globally and also agrees with the findings of Tripathi et al. (2013), who detected 6 OTUs as the most abundant archaea in 27 tropical soils.

Archaeal alpha diversity was very low in comparison with the bacterial diversity in the same black soil samples (Liu et al., 2014). Lower diversities of soil archaeal communities were also reported in other studies (Auguet et al., 2010; Singh et al., 2012; Tripathi et al., 2013). Soil pH was commonly regarded to be the major factor in the determination of the soil bacterial diversity (Fierer and Jackson, 2006; Chu et al., 2010). Similarly, several studies also indicated that the alpha diversity of the archaea in tropical and temperate soils was negatively correlated with the soil pH (Tripathi et al., 2013, 2015). In contrast, the archaeal alpha diversity in the soils from the Tibetan Plateau was significantly negatively correlated with the soil moisture, total N, total organic C (p < 0.01), total dissolved N, NH4+-N and NO3--N but not with the soil pH (Shi et al., 2016). Furthermore, Singh et al. (2012) showed that the archaeal diversity in Mt. Fuji soils had a strong positive correlation with the soil NH4+-N, K and NO3--N. In this study, we found that the OTU richness of the archaea had no relationship with any of the soil properties measured, but the OTU richness within Thaumarchaeota and Euryarchaeota was slightly correlated with soil pH (Figure 2). This change in the tendency of the OTU richness of Thaumarchaeota is not consistent with the finding of Lu et al. (2017), who detected that the OTU richness of Thaumarchaeota in forest and cropping soils had no relationship with soil pH. However, the finding of the OTU richness of Euryarchaeota in this study was supported by the results of Hu et al. (2013) who found that the OTU richness of Euryarchaeota was significantly positively correlated with the soil pH in both the uplands and paddy soils (P < 0.001). Thus, the effects of environmental factors on the alpha diversity of archaea are very complex and require further investigation.

Within the most abundant phylum Thaumarchaeota, the responses of Nitososphaerales and Thaumarchaeota_norank at the order level, and group 1.1a-associated and group 1.1b at the genus level to sampling latitude, soil pH, total C, total N were the opposite (Table 1). This finding suggests the physiological and ecological characteristics of the different taxa of Thaumarchaeota varied greatly. Recently, Lu et al. (2017) observed that the relative abundance of Thaumarchaeota SCG (group 1.1b) and TG (group 1.1c) was positively and negatively correlated with soil pH, respectively. The tendency of the change of group 1.1b is similar to that in this study, while the relative abundance of group 1.1c in the black soil showed no significant relationship with soil pH (Table 1). The lack of a significant change of group 1.1c with the soil pH in black soils could be related to its very low abundance (average 0.12%), but it accounted for ~25% of the Thaumarchaeota sequences in the study of Lu et al. (2017).

The relative abundances of the seven most abundant OTUs of Thaumarchaeota responded strongly to a number of environmental variables (Table 1). Three OTUs (OTU1, 3 and 5) were significantly correlated to the sampling latitudes, suggesting that the biogeographic distribution of archaea along the latitude exist in the black soil region. Additionally, three OTUs (OTU2, 3 and 1120) are negatively and OTU6 positively correlated with the soil pH, suggesting that the different OTUs are niche-specialized for growth at various pH levels (Tripathi et al., 2013). This finding supports the result that the soil pH is an important ecological niche to determine not only the AOA distribution (Nicol et al., 2008; Gubry-Rangin et al., 2011) but also the archaeal distribution (Tripathi et al., 2013, 2015).

The biogeographic distribution of the archaeal communities at a large scale is controversial. Ma et al. (2017) investigated the biogeographic patterns of soil bacterial, fungal, and archaeal communities across four vegetation forest zones in Eastern China, and they found that geographic distance had no influence on the archaeal communities, while the bacterial and fungal communities showing a distance-decay pattern. Similarly, (Liu et al., 2016b) also indicated that the archaeal community compositions in the lake sediments of Tibetan Plateau were not correlated with the geographic distance or altitude. Zhang et al. (2018) recently found that the bacterial, archaeal and methanogenic communities in the Chinese wetlands from all different regions did not show a significant distance-decay pattern. In contrast, Shi et al. (2016) revealed that the soil archaeal communities in the Tibetan Plateau displayed a strong distance-decay. In this study, similar to the findings of bacterial communities in the same soil samples (Liu et al., 2014), we also observed that the archaeal communities were geographically distributed across the black soil zone. However, the geographic distance on the explanation of the archaeal and bacterial community distribution was difference between two studies. The geographical distance explained 19% of the variation of the archaeal communities, which was larger than that of the bacterial communities (14.75%) (Liu et al., 2014) but was comparable to the AOA communities (17.01%) (Liu et al., 2018). Since the majority of archaea (Thaumarchaeota group 1.1b and group 1.1a) obtained in this study were phylogenetically grouped to be AOA (Nicol and Schleper, 2006), the biogeographic distribution of the archaeal communities in the black soils, to some extent, was somewhat consistent with the AOA communities as we reported previously (Liu et al., 2018).

Similar to the findings of the biogeographical distribution of the soil bacterial communities (Lauber et al., 2009; Rousk et al., 2010; Griffiths et al., 2011), pH has been shown to be the dominating soil factor in driving archaeal communities in various ecosystems (Cao et al., 2012; Hu et al., 2013; Tripathi et al., 2013, 2015). In contrast, (Liu et al., 2016b) found that the sediment salinity was the sole factor in determining archaeal communities across lake sediments in the Tibetan Plateau. In comparison, Zheng et al. (2013) reported that the soil pH, sample depth and longitude played a key role in shaping archaeal distribution in the non-flooded soil habitat, while sampling depth, longitude and NH4+-N were the most important factors in the flooded soil habitat. Recently, Zhang et al. (2018) found that the biogeographic patterns of the archaeal communities along the latitudinal gradient in the Chinese wetland soils were driven by multiple factors, including the soil pH, total organic carbon, total phosphorus, mean annual temperature, annual frost days, mean annual precipitation, and direct solar radiation. Our study showed that the distribution patterns of some archaeal taxa were significantly correlated with some soil factors (Table 1), and the soil pH was the dominant factor (contributed 20.89% variation) in the determination of the archaeal community distribution across the black soils. This finding is consistent with our previous studies showing that the soil pH was the key edaphic factor in the determination of the bacterial and AOA communities in the same soil samples (Liu et al., 2014, 2018).

In conclusion, this study revealed that the abundance of archaea in the black soils occupied 3.11~7.80% of the total prokaryotes. Multiple soil factors were found to regulate the archaeal abundance in the black soils that differed from bacteria and fungi (Liu et al., 2015; Liu et al., 2016b). Among the three observed phyla, Thaumarchaeota was the dominant archaea, followed by Euryarchaeota, and Crenarchaeota was occasionally observed at a low abundance. This study also suggests that the composition of the archaeal communities (with only 185 OTUs observed) is simplified compared to the bacterial or fungal communities in the same soils. The opposite effects of the environmental factors on the dominant OTUs suggest that the physiological characteristics of the different archaeal members are diversified in the black soils. Similar to the bacterial communities, the archaeal communities across the black soil zone are also geographically distributed with the soil pH being the major determining soil factor. In addition, Thaumarchaeota group 1.1b (Nitrososphaera) was the dominant genus, suggesting that the majority of archaea play important roles in nitrogen cycling in the upland black soils.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.