This study was carried out at the Haibei Alpine Meadow Ecosystem Research Station located in the northeast of the Tibetan Plateau. The region has a typical plateau continental climate characterized by cold and long winter but warm and short summer, with annual mean air temperature of -1.7°C and annual mean precipitation of 570 mm (Zhao et al., 2006). The vegetation growing season is from May to September. The major soil type is Mat Cryic Cambisol based on Chinese Soil Taxonomic Classification System (Arthur et al., 2008).

Soil samples were collected at elevations of 3200, 3400, 3600, and 3800 m asl. Geographical distances between adjacent elevations were 1.3 km (3600–3800 m), 4.2 km (3400–3600 m), and 6.2 km (3200–3400 m). To reduce soil heterogeneity, five soil cores at the 0–10 cm and 10–20 cm depths were randomly taken from a location with a size of 1 × 1 m and thoroughly mixed to make a composite sample. Three such locations at each elevation, which were several to dozens of meters away from each other, were chosen to create biological triplicates. All soil samples were collected in August 2009 and immediately transported on ice to the laboratory. After being sieved with 2 mm mesh to remove visible vegetation roots, residues and stones, parts of samples were stored at -80°C for DNA extraction and 4°C for environmental variable measurements, respectively.

After mixing soil cores from 0 to 10 cm and 10 to 20 cm depths, we extracted soil DNA with FastDNA spin kits (MP Biomedical, Carlsbad, CA, United States) following manufacturer’s instructions and then purified the extracted DNA in a mixture of 2.0 volume of 100% ice-cold ethanol and 0.1 volume of 3 M NaOAc (pH 5.2). After overnight incubation at -20°C, DNA solutions were centrifuged at 13000 g for 30 min and then washed with 1 ml 70% ethanol. Then DNA solutions were air dried for 30 min and dissolved in nuclease-free water. We used NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, United States) to determine DNA quality and a PicoGreen method using a FLUOstar Optima (BMG Labtech, Jena, Germany) to determine DNA concentrations (Ahn et al., 1996).

DNA was labeled with Dye Cy-5 by a random priming method as previously described (Yang et al., 2017) and the labeled DNA was purified with QIA quick purification kits (Qiagen, Valencia, CA, United States). SpeedVac (ThermoSavant, Milford, MA, United States) was used to dry DNA at 45°C for 45 min. We hybridized DNA with GeoChip 4.0 on a MAUI hybridization station (BioMicro, Salt Lake City, UT, United States) at 42°C for 16 h, as described previously (Yue et al., 2015). The slides were scanned using a NimbleGen MS200 scanner (Roche, Madison, WI, United States) at 633 nm using 100% laser power and 75% photomultiplier tube gain. Signal intensities were quantified by scanned images.

Sorensen similarity was calculated for generating the EDR z-value. An equation of log(S_S) = constant-2zlog(D) was used to calculate power-law exponent z, in which D was the Euclidean distance between two elevations and S_S was the pairwise similarity in microbial compositions by using Sorensen’s index. Elevation differences for samples at the same elevation were zero, which were changed into 0.01 to enable logarithm transformation (Legendre and Legendre, 1998). Since our pairwise comparison data were not independent with each other, we determined whether the slope was significant different from zero by bootstrapping (999 times) (Efron and Tibshirani, 1994). To this end, the original dataset was randomly paired (with replacement) and the slope of EDR was calculated by the new dataset. Using one sample t-test, we determined whether the slopes of the regressions based on bootstrapping differed from the observed slope followed by BH correction for P-value (Zhou et al., 2008).

We processed raw data as follows: (i) removing spots with a signal to noise ratio of less than 2.0; (ii) removing spots detected only once in three replicates; (iii) normalizing intensity of each spot by dividing total signal intensity of a sample and multiplying by a constant; and (iv) taking natural logarithm (Wu et al., 2016b). Subsequently, the detrended correspondence analysis (DCA), a multivariate statistical method (Ayoub-Hannaa et al., 2013), was used to compare overall functional gene structures across different samples. The hierarchical clustering analysis was used to cluster different microbial functional gene structures among soil samples. The dissimilarity test of adonis, one of the permutational multivariate analysis using distance matrices (Ricotta and Burrascano, 2009), was used to test differences of functional gene structures among different elevations. Alpha diversities were presented by Shannon and Simpson indexes. The canonical correspondence analysis (CCA) was used to identify main environmental drivers in shaping microbial communities. Variance inflation factors (VIF) of less than 20 were used as a threshold to remove redundant variables before performing CCA. The variation partitioning analysis (VPA) was used to partition environmental variables selected by CCA into three groups. Simple Mantel tests were performed to determine the relationships between microbial communities and environmental variables. Correlations between relative abundances of functional genes and environmental variables were determined by Pearson correlation and then adjusted by false discovery rate (FDR) method. These analyses were carried out using functions in the vegan package (v. 2.3-1) in R-studio v. 3.2.2. Association network analysis was performed with a subset of GeoChip data since the whole dataset was beyond the computational capacity. A random matrix theory (RMT)-based algorithm was used to construct networks using the coefficients of Spearman correlations to define similarity matrix (Deng et al., 2012).

We recorded in situ soil temperatures using type-K thermocouples (Campbell Scientific, Logan, UT, United States) with a CR1000 datalogger at the depth of 10 cm, which represents 0–20 cm. We measured in situ soil moisture at the depth of 10 cm using a time domain reflectometry (Model Diviner-2000, Sentek Pty Ltd., Stepney, SA, Australia). Soil pH data were missing for our sample. Therefore, we used soils nearby our plots collected in 2012 as substitutes since soil pH in Tibet was a generally stable variable, which had decreased only by 0.5 units since 1980 (Zhou et al., 2009). We measured soil pH with a combination glass electrode in soil/water solution of volume ratio 1:5. The concentrations of ammonia and nitrate were detected by a FIAstar 5000 Analyzer (Yang et al., 2013). To generate finer measurements of nitrogen, carbon and phosphorus, we examined total nitrogen (TN), total organic carbon (TOC) and total phosphorus at both depths of 0–10 cm and 10–20 cm. We determined TN by a Vario EL III Elemental Analyzer (Elementar, Hanau, Germany), TOC by a TOC-5000A analyzer (Shimadzu Corp., Kyoto, Japan) and total phosphorus by a UV-visible spectrophotometer (Agilent 8453, Agilent Co., Santa Clara, CA, United States). Using a Agilent 4890D gas chromatograph (Agilent Co., Santa Clara, CA, United States) equipped with Flame Ionization Detector (FID) and Electron Capture Detector (ECD) (Wang and Wang, 2003), we quantified greenhouse gasses of CO₂, CH₄ and N₂O except the 3600 m site, where rat holes were too rampant to perform accurate measurements of greenhouse gasses (Yang et al., 2014). During gas sampling, static stainless steel chambers (40 cm × 40 cm × 40 cm) were manually closed for half an hour. Then gas samples were taken every 10 min by plastic syringes and sampled within 24 h (Hu et al., 2010). We determined vegetation variables of species number, aboveground vegetation biomass and vegetation diversity (Shannon index).

Soil temperature (T) and pH decreased along the altitudinal gradient (Table 1), while soil moisture, nitrate (NO₃^-) and TN measured at the 10–20 cm depth (TN_10-20cm) increased with elevations. TOC measured at the 10–20 cm depth (TOC_10-20cm), TN measured at the 0–10 cm depth (TN_{0-10 cm}) and vegetation diversity were higher at 3400 and 3600 m than those at 3200 and 3800 m. Notably, there were correlations (P < 0.050) among TOC, TN, and total phosphorous (TP) (Supplementary Table S1).

A total of 52864 genes were detected by GeoChip, ranging from 31528 to 43816 for each site (Supplementary Table S2). The results of the dissimilarity test of adonis showed that samples were grouped by elevations (Supplementary Table S3). Microbial functional diversities, represented as Shannon indices, ranged in 10.27–10.58 and increased (P < 0.050) with elevations (Supplementary Table S2). However, the evenness indices of microbial functional genes were similar across elevations (Supplementary Table S2).

A significant EDR (r = -0.80, P < 0.001) was observed for the whole functional community (Figure 1) with a z-value of 0.0093 (Supplementary Table S4). We found that almost all z-values of carbon or nitrogen cycling processes were lower than 0.01, with the highest z-value (0.0127) of anammox and the lowest z-value (0.0055) of nitrogen fixation.

Ubiquitous genes, defined as genes detected in all of four elevations, accounted for 50.41% of all genes detected by GeoChip (Supplementary Figure S1). In contrast, unique genes, defined as genes detected only in one elevation, accounted for only 20.98% of all detected genes (Supplementary Figure S1). The number of unique genes increased with elevations, which was in accordance with the increasing ratio of the unique gene number to the detected gene number at each elevation (Supplementary Table S5). The percentage of shared genes between different elevations was the highest between 3600 and 3800 m asl (65.47%), which was the shortest in distance (Supplementary Table S6). In contrast, it was the lowest between 3200 and 3800 m asl (54.57%), which was the highest in distance. A negative correlation between distance and the number of shared genes (r = -0.71, P < 0.001) (Supplementary Figure S2) was detected.

We performed simple Mantel tests to examine major factors linking to overall microbial functional gene structure. The results showed that elevation was the most important factor in affecting functional genes structure (Table 2). To discriminate whether elevational influence is attributable to environmental variables, we performed CCA to determine the linkage between microbial functional gene structure and environmental variables (Supplementary Figure S3). Seven environmental variables (T, pH, NO₃^-, TP_10-20cm, aboveground vegetation biomass, vegetation species number and vegetation diversity), were identified to have a significant influence on functional microbial structure, which was also verified by Mantel tests (P < 0.05) (Table 2). The VPA showed that these seven variables explained a total of 82.64% variations of functional gene structure (Supplementary Figure S4).

Almost all of carbon cycling genes (Supplementary Figure S5) increased in relative abundance with elevations. However, there are certain exceptions. The relative abundances of pulA genes (encoding extracellular polyurethane/esterase) (Supplementary Figure S5A), and rubisco genes (encoding ribulose-1, 5-bisphosphate carboxylase/oxygenase) were lowest at the 3800 m site (Supplementary Figure S5B).

We found that soil moisture has a strong positive correlation with mcrA genes (encoding methyl coenzyme M reductase subunit A, r = 0.75, P = 0.010) and pmoA genes (encoding particulate methane monooxygenase A, r = 0.84, P = 0.001) (Figures 2A,B). In addition, the important role of microbes on CO₂ flux was verified by Mantel tests (r = 0.58, P < 0.001) (Table 2). In-depth investigation showed that CO₂ flux negatively correlated with aclB genes (encoding ATP citrate lyase involved in reductive TCA cycle, r = -0.83, P = 0.009), CODH genes (encoding carbon monoxide dehydrogenase involved in biosynthesis of acetyl-CoA, r = -0.89, P = 0.001) and pcc genes (involved in hydroxypropionate pathways, r = -0.77, P = 0.022) (Figures 2C–E).

To identify major carbon fixation bacteria responsible for CO₂ flux, we performed a molecular ecological network analysis with carbon fixation genes and environmental variables. The resulting network showed typical topological features of scale-free (power law = 0.92), small world (average path distance = 3.49), modular (modularity = 0.64) and hierarchical (average clustering coefficient = 0.21) (Supplementary Table S7). Notably, soil CO₂ flux showed correlations with rubisco genes derived from Roseovarius, Synechococcus, and Nakamurella multipartita, and pcc genes derived from Mycobacterium and Comamonas testosteroni (Figure 3A).

Almost all of nitrogen cycling genes increased in relative abundances with elevations, except that nosZ genes (encoding nitrous oxide reductase) were similar across four elevations (Supplementary Figure S6). NO₃^- positively correlated with amoA (encoding active site polypeptide of the ammonia monooxygenase)-AOA (ammonia-oxidizing archaea; r = 0.51, P = 0.087) but negatively correlated with amoA-AOB (ammonia-oxidizing bacteria; r = -0.60, P = 0.044) (Figures 2F,G). In addition, NO₃^- positively correlated with a number of nitrogen cycling genes, such as amoA genes (r = 0.57, P = 0.051), nirK/S genes (encoding a nitrite reductase, r > 0.58, P < 0.048), napA genes (encoding periplasmic nitrate reductase, r = 0.68, P = 0.024) and gdh genes (encoding glutamate dehydrogenase, r = 0.60, P = 0.046), and ammonia (NH₄⁺) negatively correlated with nirK/S genes (r < -0.58, P < 0.046) (Supplementary Table S8). The ratio of nirK to nirS gene abundances positively correlated with NH₄⁺ (r = 0.59, P = 0.046) and negatively correlated with moisture (r = -0.46, P = 0.123) (Supplementary Table S8). N₂O flux was positively correlated with nosZ gene derived from Roseobacter denitrificans OCh 114 (r = 0.71, P = 0.034) and another two probes derived from uncultured bacteria (r > 0.69, P < 0.036) (Supplementary Table S9), and negatively correlated with norB gene derived from Moraxella catarrhalis (r = -0.88, P = 0.002) (Supplementary Table S9).

The association network with amoA genes and environmental variables showed topological features of scale-free, small world, modular and hierarchical (Supplementary Table S7). Notably, there was a positive correlation between N₂O and amoA gene derived from Pseudomonas (Figure 3B), which was abundant in Tibet (Zhang et al., 2007; Tong et al., 2008). In addition, soil pH was negatively correlated with amoA gene derived from Aurantimonas.

Induced under cold environment, cold shock proteins help maintain cell phenotypes and viabilities through preventing damage of ice crystal (Willimsky et al., 1992). In our study, increased cold shock genes of cspA and cspB (encoding cold shock proteins) along the altitudinal gradient suggested stronger capability of microorganisms in acclimation of cold environments (Supplementary Figure S7A). In addition, we found that soil temperature correlated with cspA genes (r = -0.74, P = 0.010) and cspB genes (r = -0.90, P < 0.001) (Supplementary Table S8).

Relative abundances of genes associated with oxygen limitation, radiation stress and mercury generally increased with elevations (Supplementary Figures S7B,C,E), such as arcA genes (encoding the cytosolic response regulator) and arcB genes (encoding the membrane-bound senor kinase) derived from Shewanella, Vibrio, Nitrosospira, and Oxalobacter. In contrast, genes associated with glucose limitation were similar in relative abundances (Supplementary Figure S7D). We found that TOC_{10-20 cm} was positively correlated with mer genes encoding a mercury-resist operon (r = 0.85, P = 0.001) (Figure 2H). Soil pH was negatively correlated with mer (r = -0.78, P = 0.001) and merG (r = -0.64, P = 0.021) (Supplementary Table S8).

It has been shown that soil microbial biomass, taxonomic compositions, and bacterial diversities differed along the altitudinal gradient (Faoro et al., 2010; Fu et al., 2012). Our studies provided evidence that functional gene structures of soil microbial communities also differed along the altitudinal gradient (Supplementary Table S3). We observed increased soil moisture, TN_{10-20 cm} and NO₃^- at higher elevations (Table 1), which could provide more electron acceptors and available substrates to stimulate elemental cycling (Xu et al., 2014). In addition, it has been shown that water-extractable organic carbon, microbial biomass carbon, and carbon storage increased with elevations (Pabst et al., 2013; Tsui et al., 2013). Therefore, the increased pattern of microbial diversity could mainly be attributed to the more resources at higher elevations.

The diversity patterns of microbial communities in this study differed from those in enclosed grasslands without livestock grazing using the same data analysis method (Yang et al., 2014). We found microbial functional diversity increased with elevations (Supplementary Table S2). However, microbial functional diversity showed a humpback shape in enclosed grasslands, suggesting that land uses imposed a significant impact on soil microbial communities. Environmental stresses may not contribute to differences in microbial biogeographic patterns, as genes associated with cold shock also increased with elevations in the enclosed grassland. Rather, we speculate that differences in soil nutrient inputs could lead to differences in microbial communities. Animal excreta from livestock provide available nitrogen to vegetation and microbial communities and consequently stimulate nitrogen cycling (Kohler et al., 2005; Ambus et al., 2007). Aboveground vegetation grazed by livestock could reduce litter input to soil, which ameliorates nitrogen limitation for microbial communities because vegetation C/N ratios are higher. These two reasons could collectively stimulate substrate availability to soil microbial communities in pastures.

Livestock grazing could result in more homogenous soil environments because of soil compaction induced by trampling of livestock (Zhao et al., 2007). Similarly, we found that variations of environmental variables in grazing sites, such as NH₄⁺ (1.34–2.89 mg kg^-1), soil inorganic nitrogen (SIN; 2.88–3.62 mg kg^-1) and vegetation species numbers (19.46–26.26), were smaller than variations of NH₄⁺ (1.5–3.7 mg kg^-1), SIN (2.88–4.63 mg kg^-1) and vegetation species numbers (13.67–25.00) in enclosed grasslands (Table 1). At the functional gene level, the range of Shannon diversities (10.27–10.58) was less than half of previous observations at the enclosed grassland (9.95–10.65) (Yang et al., 2014), suggesting the important role of livestock grazing on microbial diversities. However, we could not exclude the possibility that some functional genes in soils might not be represented on GeoChip, which would lead to underestimation of microbial functional diversities (Zhou et al., 2015a). Notably, GeoChip is also limited in discovering novel genes in the environment and subjected to errors derived from cross hybridization.

Field simulation studies showed that soils exposed to warmer climates decreased functional diversities (Yue et al., 2015). In addition, warming decreased microbial diversities under normal precipitation conditions (Sheik et al., 2011). It is likely that dominant species would be more advantageous while rare species would extinct, resulting in decreased diversities (Fussmann et al., 2014). Alternatively but not exclusively, warming might increase nutrient availability (Rinnan et al., 2008). Since our results show that microbial functional diversities are higher in colder environments, we predict that warming in Tibetan pastures will decrease microbial functional diversities and shift soil carbon and nitrogen cycling. However, given strong metabolic flexibility and adaptability of soil microbial communities (Bardgett and van der Putten, 2014), it remains to be an open question whether warming will exert permanent effects on functional potentials of microbial communities and other parts of the biosphere (Lau and Lennon, 2011).

Unique vegetation species number increased with elevations between 1000 and 4000 m asl in the Tibetan Plateau (Vetaas and Grytnes, 2002). Similarly, we found that numbers of unique genes, which are genes detected only at one elevation, increased with elevations (Supplementary Table S5). The co-occurrence diversification of vegetation and microbes at higher elevations can be attributed to stronger environmental selective forces (Hanson et al., 2012), which was supported by revealing seven environmental variables explaining a total of 82.64% variations of functional gene structure (Supplementary Figure S4) and strong correlation between environmental variables and functional genes (Table 2). In contrast, evolutionary or ecological drift might play a minor role.

The tight linkages between vegetation and microbial communities, as indicated by CCA (Supplementary Figure S3) and Mantel tests (Table 2), can result from complex interactions between them. Vegetation root exudates and litter input could offer diverse organic resources to influence soil microbial communities (Wardle et al., 2004; Chu et al., 2014), while microbial communities provide feedbacks to vegetation through altering nutrient availability (Lau and Lennon, 2011). In this study, aboveground vegetation biomass correlated with microbial community structures (r = 0.52, P < 0.010) (Table 2), which has been observed elsewhere (Yang et al., 2013; Yue et al., 2015). In contrast, vegetation diversity was not correlated with microbial community structures (r = 0.14, P > 0.050) (Table 2), which could be explained by recent observations that vegetation diversity could predict beta but not alpha diversity of microbial communities in grasslands (Prober et al., 2015). More aboveground vegetation biomass confers more litter and root exudate to soil, which consequently changes C/N ratio of substrates for microbial utilization and alters microbial communities.

Negative correlations between CO₂ flux and carbon fixation genes (Figures 2C–E) suggested that changes in CO₂ flux might be attributed to anaerobic acetyl-CoA pathway, reverse TCA cycle and hydroxypropionate pathways. In contrast, Calvin-Benson-Bassham cycle, as indicated by microbial rubisco genes encoding an enzyme coupling oxygenolysis of ribulose-1, 5-bisphophate and CO₂ reduction, were similar along the altitudinal gradient. However, this observation does not refute the possibility that rubisco genes play a substantial role in soil CO₂ flux. Recently, microbial rubisco gene abundance and enzyme activities were shown to be high in Tibetan grasslands (Guo et al., 2015). In addition, we have identified correlations between CO₂ flux and rubisco genes derived from Roseovarius, Synechococcus, and N. multipartita by the association network analysis. Therefore, individual functional species possessing rubisco genes might contribute to variations of CO₂ flux, despite the lack of overall changes in total abundance of rubisco genes.

Higher soil moisture strengthens substrate utilization of microbial communities (Chen et al., 2007). Furthermore, TOC, readily mineralizable carbon, water-soluble organic carbon and reducing sugar carbon increase when soil water-holding capability varies in 20–60% (Hassan et al., 2015), owing to increase in microbial accessibility to soluble nutrients, especially autotrophic carbon, ammonia and nitrate (Bell et al., 2008; Van Horn et al., 2014). Accordingly, we observed an increase of carbon degradation genes (Supplementary Figure S5A).

Soil moisture also has a strong effect on net soil CH₄ uptake rate, which is the balance between CH₄ production and oxidation (Shrestha et al., 2012). As CH₄ production and oxidation are two tightly intertwined processes (Inagaki et al., 2004), it is no surprise to note an increase in relative abundances of both CH₄ production and oxidation genes with elevations and soil moisture (Supplementary Figure S5C). Similarly, more CH₄ and available substrates induced by soil moisture can stimulate microbial methanotrophs as well as CH₄ oxidation gene abundance (Shrestha et al., 2012). Conversely, CH₄ oxidation stimulates the growth of methanogenic bacteria by supplying energy (Hoehler et al., 1994).

Many studies used nirS and (or) nirK genes as functional markers to profile diversities or structures of denitrifiers because nitrite reductase is the rate-limiting step of denitrification (Braker et al., 1998; Liu et al., 2003). The ratio of nirK to nirS gene abundance negatively and positively correlated with soil moisture (r = -0.46, P = 0.123) and NH₄⁺ (r = 0.59, P = 0.046), respectively (Supplementary Table S8), which was consistent with a recent finding that water logging decreased the nirK/nirS ratio (Yoshida et al., 2009). In addition, both NO₃^- and NH₄⁺ are associated with community compositions of nirS-type denitrifiers (Zheng et al., 2015).

The correlation between amoA and NO₃^- (r = 0.57, P = 0.051) has been similarly documented (Ding et al., 2015). Intriguingly, amoA-AOB negatively correlated with NO₃^- (r = -0.60, P = 0.044) and amoA-AOA positively correlated with NO₃^- (r = 0.51, P = 0.087) (Figures 2F,G), which could be attributed to different roles of ammonia monooxygenase, physiological pathways and electron transport mechanisms in AOA and AOB (Hatzenpichler, 2012). Consistently, it was shown in groundwater that nitrate positively correlated with copy numbers of amoA-AOA and negatively correlated with those of amoA-AOB (Qin et al., 2014). The positive correlations between NO₃^- and hzo and napA genes provide further evidence to the finding that NO₃^- can regulate the anammox and dissimilarity nitrate reduction processes (Silver et al., 2001; Trimmer et al., 2005).

Livestock manure was one of heavy metal sources to soil, including mercury (Nicholson et al., 2003; Wang et al., 2016). At higher elevations, decreased turnover rate of soil organic carbon could accelerate mercury accumulation and cold temperature could decrease mercury evasion, causing the increasing abundance of mer genes (Zhou et al., 2015b). Consistently, mer gene positively correlated with TOC_{10-20 cm} (r = 0.85, P = 0.001) (Figure 2H), which might be attributed to mercury affinity to organic matter (Zhou et al., 2015b). In addition, the negative correlations between pH and mer (r = -0.78, P = 0.001) and merG (r = -0.64, P = 0.021) genes indicated that soil acidification constrained Hg volatilization (Supplementary Table S8).