Contrasting Biogeographic Patterns of Bacterial and Archaeal Diversity in the Top- and Subsoils of Temperate Grasslands.

Biogeographic patterns and drivers of soil microbial diversity have been extensively studied in the past few decades. However, most research has focused on the topsoil, while the subsoil is assumed to have microbial diversity patterns similar to those of the topsoil. Here we compared patterns and drivers of microbial alpha and beta diversity in and between topsoils (0 to 10 cm) and subsoils (30 to 50 cm) of temperate grasslands in Inner Mongolia of China, covering an ∼1,500-km transect along an aridity gradient. Counter to the conventional assumption, we find contrasting biogeographic patterns of diversity and influencing factors for different bacterial and archaeal groups and between depths. While bacterial diversity remains constant or increases with increasing aridity in topsoil and decreases in subsoil, archaeal diversity decreases in topsoil and remains constant in subsoil. Microbial diversity in the topsoil is most strongly influenced by aboveground vegetation and contemporary climate but is most strongly influenced by the factor historical temperature anomaly since the Last Glacial Maximum (LGM) and by soil pH in the subsoil. Moreover, the biogeographic patterns of topsoil-subsoil community dissimilarities vary for different microbial groups and are overall most strongly influenced by soil fertility differences between depths for bacteria and by contemporary climate for archaea. These findings suggest that diversity patterns observed in the topsoil may not be readily applied to the subsoil horizons. For the subsoil in particular, historical climate plays a vital role in the spatial variation of bacterial diversity. Overall, our study provides novel information for understanding and predicting soil microbial diversity patterns at depth.IMPORTANCE Exploring the biogeographic patterns of soil microbial diversity is critical for understanding mechanisms underlying the response of soil processes to climate change. Using top- and subsoils from an ∼1,500-km temperate grassland transect, we find divergent patterns of microbial diversity and its determinants in the topsoil versus the subsoil. Furthermore, we find important and direct legacy effects of historical climate change on the microbial diversity of subsoil yet indirect effects on topsoil. Our findings challenge the conventional assumption of similar geographic patterns of soil microbial diversity along soil profiles and help to improve our understanding of how soil microbial communities may respond to future climate change in different regions with various climate histories.

INTRODUCTION

Soil microbes play an indispensable role in soil formation and biogeochemical cycles and hence provide key ecosystem services, including the mediation of greenhouse gas emissions and climate change (1, 2). Exploring the biogeographic patterns of soil microbial diversity is critical for understanding mechanisms underlying the responses of soil processes to climate change. Subsoil (i.e., soils residing >20 cm below ground) contains more than half of soil organic carbon (OC) globally (3). Recent experimental studies have indicated that subsoil may show different responses to global climate changes than topsoil (4, 5) due to distinct soil environments, microbial assemblages, and their functional responses to climate changes (6). Yet, studies on soil microbial diversity have focused mostly on the topsoil, while the biogeographic patterns of microbial diversity in the subsoil on a large scale remain elusive. With a soil physical environment and microbial communities that are unique compared to those of the topsoil (3), the subsoil may show divergent patterns of microbial diversity from the topsoil. Microbial community composition differences between top- and subsoils may explain their differences in soil processes and responses to global changes (6). Hence, comparing biogeographic patterns and drivers of soil microbial diversity at different depths is important to improve our understanding of soil processes in a changing world.

Soil microbial diversity is influenced by a wide array of variables, including edaphic properties (e.g., soil pH and nutrients) (7–10), vegetation (11, 12), contemporary climate (13–15), and historical climate change (16–18), etc. These variables may have differential controls on microbial diversity in the subsoil than in the topsoil due to the varied ranges and different orders of importance of these factors. For instance, linkages between vegetation and soil microbes can be directly mediated by plant species-specific symbioses or rhizodeposition (19). Given the predominant distribution of plant roots in the topsoil, the distribution and diversity of subsoil microbes may be less affected by vegetation than those of the topsoil counterparts. Similarly, contemporary climate, including precipitation (or aridity) and temperature, has been shown to have a considerable effect on the topsoil microbial diversity (13–15), by restricting microbial access to soil nutrients or moisture (19, 20) and/or accelerating metabolic rates and biochemical processes (14). Such effects, however, may be dampened at depth (21) because microbial communities have a much longer turnover time in the subsoil (22, 23) and are considered to be less affected by contemporary climates.

Historical climate change since the Last Glacial Maximum (LGM; i.e., the most recent glaciation, ca. 21,000 to 18,000 years before present) is found to be a better predictor of species richness than contemporary climate for vertebrates (24) and plants (25, 26) in Europe and North America. A recent study also suggests that climate change since the LGM may influence soil bacterial richness and composition (16). Due to the long residence time of both soil organic matter and microbial communities at depth (21, 27, 28), microbial diversity in the subsoil may be more strongly influenced by historical climate change than that in the topsoil. However, in comparison to edaphic and contemporary climatic factors, the effect of historical climate change on soil microbial diversity patterns remains poorly understood.

In addition to varied environmental influences, the diversity pattern along soil depth may vary among different microbial clades (29). Declining carbon substrate availability with soil depth leads to an oligotrophic environment at depth, which may restrict bacterial activity and promote subsurface-dwelling groups capable of utilizing recalcitrant carbon sources (30, 31). Previous studies have shown that soil bacterial diversity is typically highest in the topsoil and decreases with soil depth (29, 32). However, the diversity or abundance of different bacterial phyla may decrease (31, 32), increase (33), or remain consistent (29) along soil profiles. In comparison with soil bacterial diversity patterns, soil archaeal diversity patterns are much less explored at depth. Some studies have revealed that the relative abundance of archaea or the ratio of archaea to bacteria tends to increase with soil depth (34), while other studies have assumed that archaeal diversity decreases or remains constant along soil profiles (35–37). Hence, diversity and composition variations of different microbial groups also need to be compared to understand the mechanisms driving microbial diversity patterns at different soil depths.

Here, using amplicon-based sequencing of 16S rRNA genes, we compare the biogeographic patterns of bacterial and archaeal diversity in and between the topsoil (0 to 10 cm) and subsoil (30 to 50 cm) along a temperate grassland transect in Inner Mongolia of China. As an integral part of the Eurasian steppe, this transect spans arid to mesic ecosystems along an aridity gradient from Northeast China toward the West, covering a broad range of climates, soil physicochemical conditions, and plant species richness. Coupled with a comprehensive data set of edaphic, vegetation, and climatic (both contemporary climate and historical climate change) variables, we evaluate the relative importance of different environmental factors driving microbial diversity at different soil depths and the microbial community dissimilarity between the top- and subsoils. This study aims to test the following three hypotheses. (i) The biogeographic patterns of diversity vary between bacteria and archaea and among different groups. (ii) Microbial diversity patterns in the subsoil do not entirely mimic those in the topsoil, and the microbial community dissimilarity between the top- and subsoils varies with environmental gradients. (iii) Microbial diversity is strongly influenced by contemporary climate and vegetation in the topsoil and by historical climate change in the subsoil.

RESULTS

Geographic variations in soil bacterial and archaeal diversity.

A survey of high-throughput amplicon sequencing for the 16S rRNA was performed to cover a large portion of bacterial and archaeal domains. According to the rarefaction results (see Fig. S1 in the supplemental material), curves of soil bacterial and archaeal communities almost reached an asymptote, suggesting that the sequencing depth was appropriate for surveying most soil bacteria and archaea. After quality filtering, denoising, and removal of potential chimeras, a total of 3,531,946 and 4,086,723 high-quality sequences (grouping into 23,458 and 3,152 operational taxonomic units [OTUs] at 97% sequence similarity per sample) were obtained for bacteria and archaea, respectively.

Rarefaction curves of all samples for bacteria and archaea generated from resampled OTUs. OTUs were defined at 97% sequence similarity, and curves of soil bacterial and archaeal communities almost reached an asymptote, suggesting that the sequencing depths were appropriate for surveying most soil bacteria and archaea. MS, meadow steppe, 8 sites; TS, typical steppe, 15 sites; DS, desert steppe, 9 sites; T, topsoil, 0 to 10 cm; S, subsoil, 30 to 50 cm. Download FIG S1, TIF file, 1.8 MB.

Soil bacteria were dominated by 3 classes (Alphaproteobacteria, Betaproteobacteria, Gammaproteobacteria) and 10 phyla, including Actinobacteria, Acidobacteria, Firmicutes, Bacteroidetes, Planctomycetes, Verrucomicrobia, Gemmatimonadetes, Nitrospirae, Chloroflexi, and Armatimonadetes (Fig. S2). Among them, Actinobacteria, Alphaproteobacteria, Acidobacteria, Chloroflexi, Nitrospirae, and Verrucomicrobia are predominantly oligotrophic while Bacteroidetes, Gemmatimonadetes, Betaproteobacteria, and Firmicutes are copiotrophic (38–40). Soil archaea were dominated by three phyla, including Crenarchaeota, Parvarchaeota, and Euryarchaeota (Fig. S2). Among them, Crenarchaeota function as ammonia-oxidizing archaea (AOA) (41). Parvarchaeota are known as acidophilic (42), while Euryarchaeota function as methanogens and denitrifiers (2, 41).

Relative abundances of the dominant phyla of soil bacteria and archaea in top- and subsoils. Download FIG S2, TIF file, 2.6 MB.

The OTU richness and phylogenetic diversity (PD) of bacteria exhibited no longitudinal trends in the topsoil from southwest to northeast, while the Shannon-Wiener diversity (here called Shannon diversity) significantly decreased (r = −0.39, P = 0.026) (Fig. 1 and Fig. S3). In contrast to topsoil bacteria, the OTU richness, PD, and Shannon diversity of topsoil archaea increased along the same gradient (r = 0.44, P = 0.012, for OTU richness; r = 0.41, P = 0.022, for PD; r = 0.58, P = 0.001, for Shannon diversity) (Fig. 1 and Fig. S3). For different clades, the OTU richness, PD, and Shannon diversity decreased from southwest to northeast in the topsoil for most oligotrophic bacterial clades (P < 0.05, except for the OTU richness and PD of Acidobacteria and Verrucomicrobia) and increased for most copiotrophic bacterial clades (P < 0.05, except for the OTU richness and PD of Gemmatimonadetes) (Fig. 2). The OTU richness, PD, and Shannon diversity of rare and unclassified archaeal clades and the Shannon diversity of Parvarchaeota increased (P < 0.05) along the same geographic direction, while the Shannon diversity of Crenarchaeota decreased (P < 0.05) (Fig. 2).

FIG 1

Sampling sites and geographic variation in soil bacterial and archaeal alpha diversity and community dissimilarity. (a) Spatial distribution of sampling sites across the temperature grasslands in Inner Mongolia; (b to g) changes in bacterial and archaeal OTU richness (b and c), phylogenetic diversity (PD) (d and e), and Shannon diversity (f and g) in topsoil and subsoil with longitude; (h and i) changes in the Bray-Curtis (h) and weighted UniFrac dissimilarities (i) between topsoil and subsoil with longitude. Land cover classification is based on the Global Land Cover Characteristics Database v2.0 (https://edcftp.cr.usgs.gov/project/glcc/globdoc2_0.html). Solid lines indicate significant linear regressions (P < 0.05).

FIG 2

Geographic variation in the alpha diversity (a to r) and community dissimilarity (s to x) of different soil bacterial and archaeal functional groups. (a to i) Topsoil; (j to r) subsoil. For bacteria, two classes (Alphaproteobacteria and Betaproteobacteria) and eight dominant phyla (the rest of those listed) were categorized as oligotrophic and copiotrophic clades, respectively. Archaea included Crenarchaeota (frequently functioning in ammonia-oxidizing processes), Euryarchaeota (frequently functioning in methane generation processes), Parvarchaeota, and the rare and unclassified clades. Solid lines indicate significant linear regressions (P < 0.05; n = 32).

Geographic variation in soil bacterial and archaeal alpha diversity in topsoil and subsoils, and community dissimilarity between top- and subsoils in the study area. Download FIG S3, TIF file, 2.7 MB.

In the subsoil, the OTU richness (r = 0.46, P = 0.008), PD (r = 0.39, P = 0.030), and Shannon diversity (r = 0.44, P = 0.013) of bacteria displayed overall opposite geographic trends from those in the topsoil, increasing from southwest to northeast (Fig. 1 and Fig. S3). Similarly, the three diversity indices of most bacterial (both oligotrophic or copiotrophic) clades also increased in the same geographic direction (P < 0.05) (Fig. 2), while only one oligotrophic clade (Verrucomicrobia) showed a decline in the Shannon diversity (P < 0.05). For archaea, only the Shannon diversity of Crenarchaeota decreased from southwest to northeast (P < 0.05); the other archaeal clades showed no trends in diversity (Fig. 2).

The topsoil-subsoil Bray-Curtis dissimilarity was overall greater than the weighted UniFrac dissimilarity for both bacterial and archaeal communities, whereas the two indices showed consistent geographic patterns for bacteria and archaea (Fig. 1). Particularly, the topsoil-subsoil dissimilarity for bacterial communities did not show clear trends from southwest to northeast, while those for archaeal communities increased (r = 0.66, P < 0.0001, for Bray-Curtis; r = 0.55, P = 0.001, for weighted UniFrac) (Fig. 1 and Fig. S3). For different bacterial clades, two oligotrophic clades (Alphaproteobacteria, Verrucomicrobia) and one copiotrophic clade (Bacteroidetes) showed decreasing trends in the topsoil-subsoil community dissimilarities along the same gradient, whereas one oligotrophic (Nitrospirae) and one copiotrophic (Gemmatimonadetes) clade showed an increasing trend in community dissimilarity (Fig. 2). For archaeal clades, only the rare and unclassified archaea showed an increasing trend in the topsoil-subsoil Bray-Curtis dissimilarity from southwest to northeast (Fig. 2). Other bacterial and archaeal clades showed no significant trends (Fig. 2).

Explanatory variables for soil bacterial and archaeal diversity variations.

The OTU richness, PD, and Shannon diversity of soil bacteria and archaea in the topsoil and subsoil displayed opposite correlations with the majority of the six environmental variables, including historical temperature anomaly, contemporary climate, vegetation, soil fertility, soil pH, and soil mineral content (Table 1). Interestingly, the oligotrophic and copiotrophic bacterial clades in the topsoil showed opposite correlations with the majority of environmental variables in terms of OTU richness, PD, and Shannon diversity. The same also applied to Crenarchaeota and the remaining archaeal phyla (Data Set S1). In the subsoil, most oligotrophic and copiotrophic bacterial clades showed relatively similar correlations with environmental variables, while archaeal clades showed no consistent correlation with environmental variables (Data Set S1).

TABLE 1

Pearson correlations between soil microbial alpha diversity in the top- and subsoils and environmental variables

Variable	Pearson correlation
	Bacteria						Archaea
	OTUrichness		PD		Shannondiversity		OTU richness		PD		Shannondiversity
	Top	Sub	Top	Sub	Top	Sub	Top	Sub	Top	Sub	Top	Sub
Historical temperature anomaly	0.18	0.44*	0.24	0.37*	−0.28	0.41*	0.35*	−0.08	0.24	−0.09	0.53**	−0.07
Contemporary climate	−0.09	0.30	−0.07	0.23	−0.53**	0.24	0.51**	0.02	0.37*	−0.04	0.55**	0.01
Vegetation	−0.09	0.31	−0.08	0.23	−0.56***	0.28	0.45**	0.00	0.32	−0.07	0.61***	0.06
Soil fertility	−0.18	0.03	−0.13	0.01	−0.47**	0.02	0.39*	0.20	0.42*	0.18	0.19	0.05
Soil pH	0.06	0.02	0.04	0.06	0.41*	0.10	−0.32	0.31	−0.19	0.38*	−0.56***	−0.02
Soil mineral	−0.25	0.28	−0.19	0.27	−0.47**	0.24	0.18	0.34	0.25	0.37*	0.08	0.06

The one, two, and three asterisks after values in bold indicate significant correlations at a P level of <0.05, <0.01, and <0.001, respectively. Top, topsoil; Sub, subsoil.

Pearson correlations between soil microbial diversity and environmental variables in and between top- and subsoils. Download Table S1, DOCX file, 0.03 MB.

Pearson correlations between soil microbial diversity and environmental variables for different phyla in topsoil and subsoil. Download Data Set S1, XLSX file, 0.1 MB.

Using hierarchical partitioning, we found that the environmental variables explained 25.5% and 29.3% of the variation in the topsoil bacterial OTU richness and PD, with historical temperature anomaly having the largest independent effects, explaining 13.4% and 18.3% of the variation (P < 0.05) (Fig. 3). The same environmental variables explained 27.4% and 19.9% of the variation in the topsoil archaeal OTU richness and PD, with contemporary climate and soil fertility having the largest independent effects, respectively (P < 0.05) (Fig. 3), and showing positive relationships with archaeal OTU richness and PD (P < 0.05) (Fig. 4). In comparison with OTU richness and PD, more variation in the topsoil bacterial and archaeal Shannon diversity indices were explained by environmental variables (45.6% for bacteria and 53.4% for archaea) (Fig. 3). Among these environmental variables, vegetation had the largest independent effect, explaining 11.2% and 15.5% of the variation in bacterial and archaeal Shannon diversity indices, respectively (P < 0.05) (Fig. 3), and showing a negative and positive relationship with bacterial and archaeal Shannon diversity indices, respectively (P < 0.05) (Fig. 4). Soil mineral content and soil pH also had similar or secondary negative independent effects on the topsoil bacterial and archaeal Shannon diversity indices, respectively, explaining 11.2% and 11.0% of the variation (P < 0.05) (Fig. 3). Besides, historical temperature anomaly and contemporary climate showed negative independent effects on the topsoil archaeal Shannon diversity, explaining 10.3% and 10.1% of the variation (P < 0.05) (Fig. 3). The other environmental variables had no significant independent effects on bacterial or archaeal OTU richness, PD, or Shannon diversity (Fig. 3).

FIG 3

Relative importance of different environment variables for alpha diversity and community dissimilarity of soil bacterial and archaeal communities. (a) Bacteria; (b) archaea. The bacterial and archaeal alpha diversity in top- and subsoils was represented by OTU richness, phylogenetic diversity (PD), and Shannon diversity, while the community dissimilarity between top- and subsoils was represented by the Bray-Curtis and weighted UniFrac dissimilarities. The relative importance of different environment variables was calculated as their independent effects using hierarchical partitioning (see Table S2 in the supplemental material). The asterisks indicate significant independent effects (P < 0.05; n = 32).

FIG 4

Changes in alpha diversity and community dissimilarity of soil bacterial and archaeal communities with dominant environmental factors. The purple, salmon, and blue points represent meadow steppe (MS), typical steppe (TS), and desert steppe (DS). Solid lines indicate significant linear regressions (P < 0.05; n = 32).

First principal component (PC 1) extracted from the principal component analysis (PCA) for the five environmental groups and results of the Kaise-Meyer-Olkin (KMO) test and the Bartlett test of sphericity (BS) for variables used for PCA in the paper. The plus (+) and minus (–) signs in parentheses refer to the positive and negative correlations between the PC 1 and the individual variables. Download Table S2, DOCX file, 0.02 MB.

In the subsoil, the environmental variables explained 36.4%, 28.7%, and 41.2% of the variation in bacterial OTU richness, PD, and Shannon diversity, respectively, with historical temperature anomaly consistently having the highest positive effect (13.3%, 10.1%, and 13.1%, P < 0.05), whereas other environmental variables had no significant independent effects (Fig. 3 and 4). The same environmental variables explained 30.5%, 32.7%, and 3.8% of the variation in subsoil archaeal OTU richness, PD, and Shannon diversity, respectively, with soil pH having the highest positive effects on archaeal OTU richness (13.6%, P < 0.05) and PD (16.1%, P < 0.05) and no environmental variables exerting significant effects on archaeal Shannon diversity (Fig. 3 and 4).

The biogeographic patterns in the topsoil-subsoil Bray-Curtis and weighted UniFrac dissimilarities of bacterial and archaeal communities were dominantly influenced by environmental variables different from those influencing diversity in the top- and subsoils (Table 2). Overall, the environmental variables explained 26.1% and 36.5% of the variation in Bray-Curtis and weighted UniFrac dissimilarities for bacterial communities and 68.5% and 71.8% for archaeal communities, respectively (Fig. 3). Among them, soil fertility had the highest negative independent effects on bacterial Bray-Curtis (17.9%, P < 0.05) and weighted UniFrac (23.9%, P < 0.05) dissimilarities, while other environmental variables had no significant independent effects (Fig. 3 and 4). Contemporary climate had the highest negative independent effect on archaeal Bray-Curtis (27.9%, P < 0.05) and weighted UniFrac (27.9%, P < 0.05) dissimilarities, followed by vegetation (22.4%, P < 0.05, for archaeal Bray-Curtis; 26.4%, P < 0.05, for archaeal weighted UniFrac) and historical temperature anomaly (11.0%, P < 0.05, for archaeal Bray-Curtis) (Fig. 3 and 4).

TABLE 2

Pearson correlations of community dissimilarity between the top- and subsoils with environmental variables

Variable	Pearson correlation
	Bacteria		Archaea
	Bray-Curtisdissimilarity	Weighted UniFracdissimilarity	Bray-Curtisdissimilarity	Weighted UniFracdissimilarity
Historical temperature anomaly	−0.03	−0.22	0.59***	0.42*
Contemporary climate	−0.07	−0.31	0.80***	0.78***
Vegetation	−0.03	−0.29	0.77***	0.78***
Soil fertility difference	−0.36*	−0.53*	0.39*	0.44*
Soil pH difference	−0.16	−0.13	−0.01	0.08
Soil mineral difference	0.02	−0.06	0.17	0.15

The one, two, and three asterisks after values in bold indicate significant correlations at a P level of <0.05, <0.01, and <0.001, respectively.

Cascading environmental effects on bacterial and archaeal diversity.

Building on the above-mentioned correlation analysis, we used structural equation models to delineate the causal effects of environmental variables on soil bacterial and archaeal diversity. The structure equation models (SEMs) are developed from a priori models (Fig. 5, solid and dotted lines) based on knowledge (13, 16, 43), with potential flows of causality from all categories of environmental variables to the dependent soil bacterial and archaeal diversity. The validated SEMs yield good model fits, indicated by nonsignificant χ2 tests (P > 0.05), high comparative fit indices (CFI > 0.95), low root mean square errors of approximation (RMSEA < 0.05 and/or RMSEA < 0.08), and low Akaike information criteria (AIC) (44).

FIG 5

Structural equation models disentangling major pathways of environmental influences on soil bacterial and archaeal alpha diversity and community dissimilarity. The bacterial and archaeal alpha diversity in topsoil (left column) and subsoil (right column) was represented by OTU richness (a and b), phylogenetic diversity (PD) (c and d), and Shannon diversity (e and f). The community dissimilarity between top- and subsoils was represented by Bray-Curtis (g) and weighted UniFrac (h) dissimilarities. Black and red arrows indicate positive and negative effects (P < 0.05), respectively, and their width is proportional to their standardized path coefficients (numbers on the arrows). Gray dotted and solid arrows indicate insignificant pathways included in the a priori and final models, respectively. Black double-sided arrows indicate Pearson correlations. R2 indicates the variance of bacterial and archaeal diversity explained by the models.

In the topsoil, the constructed SEMs explained 3.2% and 25.9% of the variation in bacterial and archaeal OTU richness (Fig. 5). Among the examined environmental variables, no environmental variables had significantly direct effects on bacterial OTU richness (Fig. 5). In contrast, contemporary climate had a direct and positive effect on archaeal OTU richness (Fig. 5). The SEMs explained 26.5% and 16.6% of the variation in bacterial and archaeal PD, with historical temperature anomaly and contemporary climate having direct effects on bacterial PD and soil fertility having a direct positive effect on archaeal PD (Fig. 5). The SEMs explained 39.7% and 47.5% of the variation in the bacterial and archaeal Shannon diversity indices, respectively (Fig. 5). Vegetation and soil mineral content had a direct negative effect on the topsoil bacterial Shannon diversity, while contemporary climate and historical temperature anomaly indirectly influenced the bacterial Shannon diversity by affecting vegetation and soil mineral content, respectively (Fig. 5). In contrast, vegetation and soil fertility had direct positive and negative effects, respectively, on the archaeal Shannon diversity (Fig. 5).

In the subsoil, the constructed SEMs explained 20.8%, 11.3%, and 37.3% of the variation in bacterial OTU richness, PD, and Shannon diversity, respectively. Among the environmental variables, historical temperature anomaly had direct positive effects on the three diversity indices. Besides, soil minerals also had a direct positive effect on the bacterial PD, and vegetation and soil pH had direct positive effects on the subsoil bacterial Shannon diversity. In contrast, the SEMs explained 22.3% and 23.7% of archaeal OTU richness and PD, respectively, with soil pH and soil mineral content having direct positive effects consistently. No acceptable model was yielded for the archaeal Shannon diversity (Fig. 5).

The SEMs explained 12.9% and 28.0% of the variation in the topsoil-subsoil Bray-Curtis and weighted dissimilarities for bacterial communities and 64.3% and 68.2% of the variation for archaeal communities, respectively (Fig. 5). Among the six environmental variables, only soil fertility difference had direct and negative effects on bacterial Bray-Curtis and weighted dissimilarities. In contrast, contemporary climate had direct effects on the archaeal Bray-Curtis dissimilarity, and vegetation and soil pH differences had direct effects on archaeal weighted dissimilarities (Fig. 5). Contemporary climate only indirectly influenced archaeal weighted dissimilarities via its effect on vegetation (Fig. 5).

DISCUSSION

While the biogeographic pattern of soil microbial diversity has long been studied in the topsoil, its variation in the subsoil compared to that in the topsoil remains largely unknown. Here, using amplicon-based sequencing of 16S rRNA genes, we show contrasting microbial diversity patterns and influencing factors in the topsoil versus the subsoil as well as changes in topsoil-subsoil microbial community dissimilarities along an aridity gradient in the temperate grasslands of Inner Mongolia. Our results reveal divergent diversity patterns among different microbial phyla and functional groups at the regional scale.

Divergent diversity patterns for various microbial groups in the topsoil.

Our study reveals divergent geographic patterns of bacterial and archaeal alpha diversity in the topsoil along the aridity gradient in the temperate grasslands of Inner Mongolia, with bacterial diversity showing none (OTU richness and PD) or decreasing (Shannon diversity) trends, and archaeal diversity (OTU richness, PD, and Shannon diversity) increasing from southwest toward northeast (Fig. 1; see Fig. S3 in the supplemental material). These results support our first hypothesis that biogeographic patterns of diversity vary between bacteria and archaea and among different lineages within these groups. Among the environmental factors, historical temperature anomaly since LGM has the greatest independent effects on the topsoil bacterial OTU richness and PD and influences bacterial diversity directly. This result is consistent with previous findings that the paleoclimate in the LGM and mid-Holocene also explains a significant proportion of the global variation in topsoil bacterial diversity (16). However, the effect of historical temperature anomaly can be detected only when its shared effects with other environmental factors are controlled for (Fig. 3 and 5), suggesting that its effect may be concealed by other factors and/or due to rapid turnover of microbial population and organic matter in the topsoil (22, 23). Unlike with bacteria, contemporary climate and soil fertility dominate the geographic patterns of topsoil archaeal OTU richness and PD, respectively (Fig. 3). This result suggests that different dimensions of archaeal diversity may be driven by different environmental factors. Among contemporary climate variables, aridity index and soil water content have the strongest positive effects on the topsoil archaeal OTU richness, mainly dominating the geographic pattern of the rare and unclassified archaeal clades (Table S1 and Data Set S1). Similarly, soil fertility influences the topsoil archaeal PD, also mainly influencing rare and unclassified archaeal clades (Data Set S1). Therefore, disentangling the biogeographic distribution of nonculturable archaea is critical for understanding the biogeography of archaeal diversity (2, 45).

Vegetation has the strongest yet opposite effects on bacterial and archaeal Shannon diversity indices (Fig. 5). Furthermore, contemporary climate, including mean annual precipitation (MAP), mean annual temperature (MAT), aridity index, and soil water content, has an indirect and minor influence via its effects on vegetation (see Fig. 5). This result is somewhat in line with the previous findings that vegetation and contemporary climate exert strong effects on the topsoil bacterial Shannon diversity in arid and semiarid temperate grasslands (15), albeit in opposite directions. Bacterial communities in the topsoil of the moisture- and N-limited Inner Mongolian grasslands (38, 46) are dominated by oligotrophic clades (such as Actinobacteria and Chloroflexi) (Fig. 2) due to their higher substrate affinities relative to those of copiotrophic clades (40). Therefore, patterns in the diversity of oligotrophic clades and their response to vegetation and contemporary climate changes may have dominated the patterns and responses of total bacterial diversity.

Among the variables representing vegetation, plant species richness has a particularly strong negative effect on the bacterial Shannon diversity (Table S1). This result stands in contrast with previous studies reporting either a positive correlation between bacterial diversity and plant diversity in the Rocky Mountains in Colorado, USA (12), or a neutral relationship in global temperate grasslands (11, 47). For clades with different functions, the Shannon diversity of most oligotrophic clades is negatively correlated with plant species richness (Data Set S1), consistent with previous studies (38, 40). In contrast, the Shannon diversity of most copiotrophic clades (such as Betaproteobacteria and Gammaproteobacteria) is positively correlated with plant species richness (Data Set S1) (47). These results together further corroborate the dominance of oligotrophic bacterial clades in relatively infertile grassland such as the studied transect.

In contrast to bacteria, net primary productivity (NPP) and aboveground biomass rather than plant species richness among the variables of vegetation dominate the patterns of the archaeal Shannon diversity (Table S1). This suggests that plant carbon inputs may lead to an increase in the topsoil archaeal Shannon diversity along the aridity gradient (9, 34). Various archaeal phyla contribute differently to the Shannon diversity of total archaea. Specifically, the increase in the archaeal Shannon diversity from southwest toward northeast is dominated by Parvarchaeota and the rare and unclassified clades (Fig. 2). The Shannon diversity of Parvarchaeota, which can degrade multiple carbon resources (e.g., starch, cellulose, and disaccharides) (42), is positively correlated with NPP and aboveground biomass (Data Set S1). By comparison, Euryarchaeota are predominantly methanogens and capable of autotrophic growth, conferring their relative independence from plant carbon inputs indicated by NPP and aboveground biomass (Data Set S1) (2, 41). In contrast, the Shannon diversity of Crenarchaeota displays a negative correlation with NPP (Data Set S1). As Crenarchaeota are capable of mixotrophic growth and assimilating carbon from oxidized inorganic compounds, i.e., carbon dioxide (CO2) or bicarbonate (HCO3−) (2), they may compete strongly with plants for N in these N-limited grasslands, thus constraining their diversification under elevated plant growth.

In addition, edaphic factors, including soil mineral content, fertility, and soil pH, also influence soil bacterial and archaeal diversity. Soil mineral content negatively influences bacterial OTU richness and Shannon diversity, mainly due to the negative response of oligotrophic bacterial clades (Data Set S1). The effects of soil mineral content may reflect the effects of historical temperature anomaly on soil formation (16). Soil fertility shows causal but opposite effects on archaeal PD and Shannon diversity in the SEMs. Its effect primarily reflects the indirect effects of contemporary climate, historical temperature anomaly, and soil mineral content on archaeal diversity (Fig. 5). Among the variables within the soil fertility group, soil total nitrogen is particularly important for archaeal diversity (Table S1), consistent with a recent study in eastern China forests (48). The negative effects of soil fertility and mineral content on archaeal Shannon diversity are attributed mainly to the negative response of Crenarchaeota (Data Set S1), which frequently function as ammonia-oxidizing archaea (41).

Contrasting patterns and drivers of microbial diversity in the subsoil versus topsoil.

Our study demonstrates contrasting geographic patterns of microbial diversity in the subsoil and topsoil of the Inner Mongolian grasslands. While bacterial diversity in the topsoil shows no clear trends or decreases from southwest toward northeast, it increases in the subsoil (Fig. 1). Similarly, archaeal diversity shows opposite patterns in the topsoil and subsoil. These results support our second hypothesis and suggest that we cannot infer the biogeographic patterns of microbial diversity in the subsoil from those in the topsoil. Furthermore, as recent studies indicate that subsoil biogeochemical processes may be strongly influenced by global climate change (4–6), it is urgent to further explore patterns as well as drivers of microbial diversity in the subsoil.

Using different statistical analyses, we also find that microbial diversity is influenced by different variables in the top- and subsoil. In contrast to the dominant control by historical temperature anomaly, vegetation, and contemporary climate on bacterial diversity in the topsoil, historical temperature anomaly is the only dominant driver of bacterial diversity in the subsoil (Fig. 3), which consistently and directly influences the three dimensions of bacterial diversity (Fig. 5). This result supports our third hypothesis. Influences of historical temperature anomaly on subsoil bacterial diversity may be associated with its direct legacy effect on the distribution of soil bacteria during the past and/or indirect effects on edaphic factors. Andam et al. (49) and Martiny (17) argued that the influence of climate conditions more than 10,000 years ago can be found in contemporary soil bacterial populations, such as Streptomyces. Previous studies also suggest that historical climate (e.g., precipitation) can affect bacterial diversity directly via its influence on enzyme sensitivity (50) or indirectly via its influence on soil properties (such as carbon stocks and quality) (16, 28, 51). More importantly, we show that historical temperature anomaly directly regulates subsoil rather than topsoil bacterial diversity (except topsoil bacterial PD) (Fig. 5). Microbial population and organic matter have a much slower turnover and longer residence time in the subsoil than in the topsoil (22, 23), potentially rendering them less susceptible to contemporary than historical climate variations. Historical temperature anomaly influences subsoil bacterial diversity mainly via its effects on the diversity of Actinobacteria, Acidobacteria, Alphaproteobacteria, Bacteroidetes, Gemmatimonadetes, Betaproteobacteria, Gammaproteobacteria, and Verrucomicrobia (Data Set S1), indicating that those clades may have experienced higher variability during the last glacial period (16, 17). Besides, plant and soil pH also have direct but weak independent effects on the subsoil bacterial Shannon diversity, which might be caused by their close relationship with other factors, such as contemporary climate.

Soil pH and soil mineral content (including contents of silt, sand, and extractable Fe) significantly explain the patterns of subsoil archaeal OTU richness and PD (Fig. 3), which is in contrast to the primary drivers of topsoil archaeal diversity (i.e., contemporary climate, soil fertility, and vegetation for the three dimensions of topsoil archaeal diversity, respectively). The fact that influences of edaphic variables (soil pH and soil minerals) outcompete those of vegetation and contemporary climate on subsoil archaeal diversity might be because root biomass has a low proportion in subsoil, as most roots usually develop in the topsoil, potentially rendering subsoil archaeal diversity less susceptible to vegetation and contemporary climate than to edaphic variables (19).

Drivers of topsoil-subsoil microbial community dissimilarity.

Given the contrasting microbial diversity patterns in the topsoil versus the subsoil, we analyze patterns of topsoil-subsoil microbial diversity dissimilarity in terms of taxonomic beta diversity (Bray-Curtis) and phylogenetic beta diversity (weighted UniFrac) among different microbial groups. The pairwise community dissimilarities (Bray-Curtis and weighted UniFrac) show contrasting geographic patterns for bacteria and archaea, corroborating our second hypothesis. Community dissimilarities did not show clear trends for bacteria from arid and semiarid grasslands in the southwest to mesic grasslands in the northeast but increased for archaea (Fig. 1). These results are consistent with previous studies (29, 36) and suggest that microbial community composition may significantly vary vertically across different regions and among different ecosystems (e.g., arid versus mesic grasslands). They also emphasize the importance of soil depth as an environmental gradient that structures soil microbial communities, especially soil archaeal communities.

To further reveal environmental drivers for the topsoil-subsoil microbial community dissimilarities, we explored whether topsoil-subsoil differences in soil parameters, together with other climatic and vegetation variables, contribute to the microbial community dissimilarities. Both hierarchical partitioning and SEM analyses indicate that soil fertility difference significantly influences bacterial Bray-Curtis and weighted UniFrac dissimilarities (Fig. 3 and 5). The strong negative effect of soil fertility on community dissimilarities of bacteria is driven mainly by the community dissimilarities of some oligotrophic clades (such as Acidobacteria and Alphaproteobacteria) and some copiotrophic clades (such as Gemmatimonadetes, Betaproteobacteria, and Firmicutes) (Fig. 2). These results support previous findings on the role of soil fertility in shaping soil bacterial community composition (52). Although historical temperature anomaly, contemporary climate, and vegetation have no significant effects on total bacterial community dissimilarities, they promote shifts in the community composition of some bacterial clades (such as Armatimonadetes, Nitrospirae, and Gemmatimonadetes) (Data Set S1).

In contrast to bacteria, the topsoil-subsoil archaeal community dissimilarities are most strongly driven by contemporary climate (Fig. 3). The effect of contemporary climate on archaeal community dissimilarities is achieved primarily through its influence on the dissimilarity of Parvarchaeota and the rare and unclassified clades (Data Set S1). In addition, vegetation and historical temperature anomaly play a secondary role in driving archaeal community dissimilarities, indicating that contemporary archaeal community composition is also influenced by recent plant carbon input (9) and long-term historical climate change (17). Most archaeal members are considered to be descendants of very old cell lineages (2) and thus are more easily influenced by historical temperature anomaly.

Conclusions.

Our results demonstrate contrasting biogeographic patterns of diversity between bacteria and archaea in the studied temperate grasslands, highlighting the varied responses of different microbial groups to environmental variations in the soil. More importantly, by comparing microbial diversities at different soil depths, we show that microbial diversity patterns in the subsoil do not mimic those in the topsoil. Until now, studies have focused primarily on microbial diversity patterns in the topsoil. Our results suggest that these studies may misrepresent the distributions and diversity variations of vast microbial communities at soil depths. It is therefore essential to add a new dimension (soil depth) to our understanding of soil microbial diversity variations along spatial gradients. Furthermore, historical temperature anomaly plays a more important and direct role in regulating bacterial diversity in both the topsoil and the subsoil. Finally, we should mention that microbial communities may vary between different seasons in the temperate ecosystems, experiencing significant seasonality. Seasonal variation was not analyzed in this study. The legacy effect of historical climate change on subsoil microbial diversity and the seasonal dynamics of soil microbial community need to be considered to better understand and predict the impacts of future climate change on soil microbial diversity.

MATERIALS AND METHODS

Study area and soil sampling.

Our study area spans an ∼1,500-km transect ranging from arid to mesic grasslands in Inner Mongolia (ca. 107.929°E to 119.970°E, ca. 39.154°N to 49.618°N) with varied climatic, edaphic, and vegetation conditions (see Fig. S4 and Data Set S2 in the supplemental material). This transect includes several vegetation types (desert steppe, typical steppe, and meadow steppe) with increasing mean annual precipitation (MAP) (ca. 165.0 to 411.5 mm) and decreasing mean annual temperature (MAT) (ca. 6.4°C to –2.3°C) from southwest toward northeast. The desert steppe is arid and low in plant species richness, dominated by perennial drought-adaptive species, including Stipa klemenzii and Stipa breviflora, etc. (53). The typical steppe has the highest coverage in Inner Mongolian, with intermediate levels of NPP and plant species richness, dominated by Stipa grandis, Stipa krylovii, and Artemisia frigida, etc. (46). The meadow steppe has the highest NPP and plant species richness, dominated by Stipa baicalensis and Leymus chinensis, etc. (46). Soil types along this transect include Calcisols, Kastanozems, and Calcic Chernozem from southwest toward northeast (46).

Geographic variation in all the environmental, aboveground vegetation, and soil properties in topsoil, subsoil, and topsoil-subsoil difference in the study area. Topsoil-subsoil difference (Diff) in soil properties were calculated as follows: topsoil-subsoil difference = (Xtop – XSub)/(Xtop + Xsub) × 100%, where X represents soil properties. Positive values indicate that the topsoil has higher concentrations in soil properties than the subsoil, whereas the negative values indicate the opposite patterns. The topsoil-subsoil difference in soil properties is displayed on the right y axis. All the solid lines and dotted lines represent significant (P < 0.05) and marginally significant (0.05 < P < 0.1, n = 32) linear regression. MAP, mean annual precipitation; MAT, mean annual temperature; plant AGB, plant aboveground biomass; plant NPP, plant net primary productivity. Download FIG S4, TIF file, 2.6 MB.

Information on sampling sites and raw data in the Inner Mongolian grasslands. Download Data Set S2, XLSX file, 0.02 MB.

Soil samples were collected from 32 randomly selected sites along the transect in August 2015. At each site, five subplots (1 m by 1 m) were set at the four corners and middle of a large plot (10 m by 10 m). Three subplots along the diagonal were randomly selected for each large plot. Within each subplot, three soil cores were taken by excavating soils from predetermined depths to a total of 100 cm using a 50-mm-diameter soil auger (54). Soils from the same depth and subplot were thoroughly mixed as a composite sample and divided into two portions. One portion was kept in an ice box and stored at –80°C immediately after being transported to the laboratory for DNA analysis, while the other portion was air-dried for physicochemical analyses. In this study, only the topsoil (0 to 10 cm) and subsoil (30 to 50 cm) samples were used, and three subplot replicates were thoroughly mixed to constitute a representative sample at each site. All soils were sieved through a 2-mm mesh, with visible roots removed before laboratory analysis. The aboveground biomass (AGB) of each species was harvested by clipping the entire aboveground part, dried at 75°C to a constant weight, and weighed separately for each subplot. The NPP of each site was estimated using data from the Numerical Terradynamic Simulation Group (NTSG) with a spatial resolution of 1 by 1 km (http://www.ntsg.umt.edu/project/modis/mod17.php).

Soil physicochemical analysis.

Total carbon (TC) and total nitrogen (TN) concentrations of soil samples were measured by combustion using an elemental analyzer (Vario EL III; Elementar, Hanau, Germany). Soil OC was calculated as total carbon minus inorganic carbon, which was analyzed volumetrically by reaction with hydrochloric acid, as previously described (55). Total phosphorus (TP) was extracted using perchloric acid-sulfuric acid (HClO4-H2SO4) digestion and measured by a colorimetric method with molybdenum blue (56). Soil pH was measured using a soil-to-water ratio of 1:2.5 (wt/vol). Soil texture was examined by laser diffraction using a Malvern Mastersizer 2000 (Malvern Instruments Ltd., UK) after removal of organic matter and calcium carbonates (55). Dithionite-extractable iron (Fed) and aluminum (Ald) were extracted from soil using the citrate-bicarbonate-dithionite (CBD) method (57) and subsequently determined on an inductively coupled plasma-atomic emission spectrometer (ICP-AES; ICAP6300, Thermo Scientific, USA).

DNA extraction and high-throughput amplicon sequencing.

DNA was extracted from soils using the MoBio PowerSoil DNA isolation kit (MoBio Laboratories, Carlsbad, CA, USA) according to the manufacturer’s protocol. DNA concentration was first assessed on 1% agarose gels and a NanoDrop 2000/2000C (NanoDrop, Germany) based on 260/280 and 260/230 nm absorbance ratios. According to the concentration, DNA was diluted to 1 ng μl−1 using sterile water to serve as a template solution.

For bacteria, the V4 region of the 16S rRNA gene was amplified with the forward primer 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and the reverse primer 806R (5′-GGACTACHVGGGTWTCTAAT-3′), generating ca. 253-bp fragments (58). The primers contain a pair of 6-bp error-correcting forward and reverse barcode sequences, respectively. For archaea, 16S rRNA genes were amplified with primer pair 1106F (5′-TTWAGTCAGGCAACGAGC-3′) and 1378R (5′-TGTGCAAGGAGCAGGGAC-3′) with a pair of 8-bp forward and reverse barcode sequences, generating ca. 280-bp fragments (59). The primer set 1106F/1378R mainly targeted methanogenic archaeal 16S rRNA genes but can still detect nonmethanogenic clades due to nonspecificity (60). All the barcodes were unique to every soil sample.

The PCR was performed in 30-μl reaction systems after mixing 15 μl of Phusion high-fidelity PCR master mix (New England Biolabs), 0.2 μM forward and reverse primers labeled with specific barcodes, and about 10 ng template DNA. Thermal cycling was repeated by use of the following procedure: initial denaturation at 98°C for 1 min, followed by 30 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s, and elongation at 72°C for 30 s, with a final step of 72°C for 5 min. At the termination of thermal cycling, PCR products were mixed with the same volume of 1× loading buffer (contained SYBR green) and used to conduct electrophoresis on a 2% agarose gel for detection. Samples with a bright main strip between 400 and 450 bp were chosen for further experiments. PCR products were mixed in equal density ratios, and then PCR mixture products were purified with a GeneJET gel extraction kit (Thermo Scientific). Equal molar concentrations of PCR products for each sample were pooled. Sequencing libraries were generated using an Illumina TruSeq DNA PCR-free library preparation kit (Illumina, USA) in accordance with the manufacturer’s recommendations, and index codes were added. The library quality was assessed using a Qubit 2.0 fluorometer (Thermo Scientific) and an Agilent Bioanalyzer 2100 system. Finally, the libraries were sequenced on an Illumina HiSeq 2500 platform, and paired-end reads were generated in fastq or fasta format with forward and reverse directions assigned to separate files.

Processing of sequencing data.

Raw DNA sequences generated from the Illumina HiSeq 2500 platform were processed on the Galaxy pipeline in Metagenomics for Environmental Microbiology (http://mem.rcees.ac.cn:8080/root/index) (61) at the Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences. Specifically, the raw DNA sequences assigned to samples were first cleaned by removing the barcodes and primer sequences. The paired-end reads were then merged by FLASH (version 1.0.0), a very fast and accurate analysis tool which is designed to merge paired-end reads (62). The minimum required overlap length of paired-end reads was set to at least 30 bp, and the maximum overlap length approximated 90% of read pairs. The maximum allowed ratio of number of mismatches to overlap length was set as 0.25, with the Phred Offset representing the quality values of bases set as 33 and the standard deviation set as 10% of the average fragment length. After merging the paired-end reads, the sequences were filtered with the BTRIM program with an average quality score threshold of >20 over a 5-bp window size and a minimum length of 200 bp (63). The sequences were further denoised by removing the sequences of less than 200 bp or with ambiguous bases. Finally, the sequences were trimmed to keep sequences for bacteria between 245 and 260 bp and for archaea between 272 and 288 bp, followed by exclusion of putative chimeric sequences. Therefore, we obtained a total of 3,531,946 and 4,086,723 high-quality bacterial and archaeal sequences, respectively, which were grouped into 23,458, and 3,152 OTUs for soil bacteria and archaea at 97% sequence similarity, and corresponding fasta format sequences were obtained using the UPARSE pipeline (64).

The reads of OTUs were annotated by referring to the Greengenes database (65) for taxonomic information for bacteria and archaea with a minimum 50% confidence score. Because the 505F/806R primer pair can target a small quantity of archaea due to its marginal nonspecificity, we therefore removed the OTUs that were annotated as archaea in the following analysis. In addition, the OTUs annotated as bacteria by the 1106F/1378R primer set targeting archaea were also removed.

To build a phylogenetic tree with the fasta sequences, MAFFT software (66) was first used to align the sequences, and a maximum likelihood (ML) tree was built using ExaML software (67) for soil bacteria and RAxML software (68) for soil archaea. Both ExaML and RAxML were obtained from https://cme.h-its.org/exelixis/software.html. To make the data comparable among different sites, we standardized the OTU table across all samples to 32,885 and 50,347 sequences (all were the smallest number of sequences across the sample) for bacteria and archaea per sample, respectively. All the following analyses were based on the standardized data. By random sampling and generation of rarefaction curves, we found that the rarefaction curves for all samples for soil bacterial and archaeal OTUs leveled off at the current sequencing depth.

Climate data.

To evaluate the effect of contemporary climate on soil microbial diversity, we used mean annual precipitation (MAP, mm), mean annual temperature (MAT, °C), aridity index, and soil water content (SWC, mm month−1) from 1950 to 2000. These variables have been shown to be the dominant factors of both aboveground and belowground communities in the Inner Mongolian grassland in previous studies (15). The MAP and MAT data, with a spatial resolution of 30 arc seconds, were obtained from the WorldClim website (http://worldclim.org/version2) (69). The aridity index is calculated as the ratio of MAP to potential evapotranspiration (PET). The data on PET and soil water content with a spatial resolution of 30 arc seconds were obtained from the CGIAR-CSI Global PET database (www.cgiar-csi.org/data/global-aridity-and-pet-database) and soil water balance database (https://cgiarcsi.community/data/global-high-resolution-soil-water-balance), respectively (70).

To calculate the climate data of a site with a given longitude, latitude, and altitude, we took the following steps. First, the grid cells of a data layer within 100 km from the site were extracted. Second, the longitude and latitude of the centroids of these grid cells were calculated, and their altitudes were extracted from the GTOPO30 digital elevational model with a resolution of 1 by 1 km (http://eros.usgs.gov/#/Find_Data/Products_and_Data_Available/gtopo30_info) using their centroid coordinates. Third, the following model was established for each variable separately, using the extracted climate data, longitude, latitude, and altitude of these grid cells:where a, b, c, and d are regression coefficients, MMT is mean monthly temperature, and MMP is mean monthly precipitation. Fourth, the value of each variable at the focal plot were calculated separately by inputting the longitude, latitude, and altitude into the corresponding model.

To evaluate the effect of historical climate change on soil microbial diversity, we calculated the anomaly of mean annual temperature (T anomaly) as contemporary mean annual temperature minus that at the Last Glacial Maximum based on MIROC (Model for Interdisciplinary Research on Climate) (71).

Statistical analysis.

OTU (operational taxonomic unit) richness, phylogenetic diversity (PD), and Shannon-Wiener diversity (here called Shannon diversity) were used to estimate the alpha diversity in topsoil and subsoil microbial communities with the vegan (version 2.4-5) package in R (version 3.4.3) (72). Bray-Curtis and weighted UniFrac dissimilarities were used to estimate the taxonomic and phylogenetic beta dissimilarities between paired topsoil and subsoil microbial communities, respectively, using vegan and phyloseq (version 1.22.3) packages in R (73).

To explore the drivers of microbial alpha and beta diversity variations, 19 variables were compiled or measured, including historical temperature anomaly since LGM, MAP, MAT, aridity index, soil water content, plant aboveground biomass, plant species richness, NPP, soil total nitrogen, soil total carbon, soil organic carbon, soil total phosphorus, soil pH, soil-extractable Ca, soil-extractable Mg, soil-extractable Fe, soil-extractable Al, soil clay, soil silt, and soil sand (Fig. S4, Table S1, and Data Set S2; see details in Text S1). To avoid collinearity between variables in the following regression analysis, we classified all parameters into six groups based on their ecological implications: (i) soil fertility (including soil total nitrogen, soil total carbon, soil organic carbon, and soil total phosphorus); (ii) soil pH; (iii) soil mineral content (including soil silt, soil sand, soil-extractable Fe, and soil-extractable Al); (iv) vegetation (including plant aboveground biomass, plant species richness, and NPP); (v) contemporary climate (including MAP, MAT, aridity index, and soil water content), and (vi) historical temperature anomaly. Differences in soil properties between the topsoil and the subsoil were calculated as follows:where X represents the soil property, top represents topsoil, and sub represents subsoil.

Supplementary results: changes in climate and soil properties along the transect. Download Text S1, DOCX file, 0.02 MB.

Principal-component analysis (PCA) was conducted for each group encompassing more than one variable, and the first principal component (PC 1) was extracted to represent each variable group. These components explained 62.2% to 92.6% of the variations in the original variables (Table S2). The feasibility of using PCA was checked using the Kaise-Meyer-Olkin (KMO) test and the Bartlett test of sphericity (BS) (Table S2), which indicates that PCA is appropriate to use for our data (43).

Relationships of microbial diversity and topsoil-subsoil community dissimilarity with environmental variables were assessed by a simple Pearson correlation using the R package Hmisc (74). To further compare the independent effects of different environmental factors, we conducted hierarchical partitioning using the R package hier.par (75). The relative independent effects refer to their independent effects in the total variations.

Structure equation models (SEMs) were used to evaluate the direct and indirect effects of environmental factors on microbial diversity and topsoil-subsoil community dissimilarity (13). The SEMs were fitted by maximum likelihood estimation using AMOS 17 (13). For the categories of environmental variables, the PC 1 of the four variable groups and two individual variables (historical temperature anomaly and soil pH, which were standardized) were used as predictors. A priori models were evaluated and optimized by stepwise exclusion of variables with nonsignificant regression weights and stepwise inclusion of additional correlations based on modification indices and goodness of fit for the initial model (43). Due to our relatively small data set with a nonnormal distribution, the models were modified with the Satorra-Bentler correlation to improve the chi-square approximation of goodness-of-fit test statistics and confirmed using the Bollen-Stine bootstrap test (43). Models were considered to have a good fit when the bootstrap P value was within 0.1 to 1.0. Since there is no single universally accepted test of overall goodness of fit for SEMs, we also used the χ2 test, the root mean square error of approximation (RMSEA), the CFI, and the AIC as criteria to test the goodness of the model fit (13). The model has a good fit when the χ2 and AIC are low, the CFI is high (CFI > 0.95), and the RMSEA is near 0 (RMSEA values of ≤0.05 can be considered a good fit; values between 0.05 and 0.08 can be considered an adequate fit) (44). We checked the bivariate relationships between all variables to ensure that a linear model was appropriate (Fig. S5).

Correlation matrix of environmental variables used in the SEMs for microbial diversity in topsoil (a) and subsoil (b) and community dissimilarity between topsoil and subsoil (c). Numbers in the lower left triangle indicate Pearson correlation coefficients (r) between the corresponding variables (P < 0.05), proportional to the size of the colored dots in the upper right triangle. Empty cells indicate nonsignificant correlations (P > 0.05). Download FIG S5, TIF file, 2.7 MB.

Data availability.

HiSeq 2500 sequencing data have been deposited in the public National Center for Biotechnology Information (NCBI) database under BioProject accession number PRJNA557316.