Long-term warming in a Mediterranean-type grassland affects soil bacterial functional potential but not bacterial taxonomic composition

Climate warming is known to impact ecosystem composition and functioning. However, it remains largely unclear how soil microbial communities respond to long-term, moderate warming. In this study, we used Illumina sequencing and microarrays (GeoChip 5.0) to analyze taxonomic and functional gene compositions of the soil microbial community after 14 years of warming (at 0.8–1.0 °C for 10 years and then 1.5–2.0 °C for 4 years) in a Californian grassland. Long-term warming had no detectable effect on the taxonomic composition of soil bacterial community, nor on any plant or abiotic soil variables. In contrast, functional gene compositions differed between warming and control for bacterial, archaeal, and fungal communities. Functional genes associated with labile carbon (C) degradation increased in relative abundance in the warming treatment, whereas those associated with recalcitrant C degradation decreased. A number of functional genes associated with nitrogen (N) cycling (e.g., denitrifying genes encoding nitrate-, nitrite-, and nitrous oxidereductases) decreased, whereas nifH gene encoding nitrogenase increased in the warming treatment. These results suggest that microbial functional potentials are more sensitive to long-term moderate warming than the taxonomic composition of microbial community.


INTRODUCTION
Accumulating concentrations of greenhouse gases in the atmosphere are increasing global surface temperature, which is expected to persist in the coming decades 1 . Warming often stimulates plant photosynthesis, biomass, and growth [2][3][4][5] , and alters plant community composition 6,7 . Soil C and N cycling are also altered by warming [8][9][10][11] . In comparison, there is less information on responses of the soil taxonomic and functional genes to warming. More precisely, many studies have analyzed the effect of brief heat shocks (often applied over a few minutes or hours) 12,13 ; or warm climate spells (typically applied over a few days or weeks) 14,15 ; on soil microbial communities. But studies analyzing the effect of moderate and progressive warming on a long term (over years to decades), as currently experienced due to global warming, are very scarce.
Soil microbial communities are vastly diverse, and are key drivers of soil biogeochemical cycling 16,17 . Warming can increase soil microbial activities, change microbial biomass, and community composition 15,[18][19][20] , but the magnitude of warming effects vary by climate regions and ecosystem type. For example, a metaanalysis of 64 studies showed that warming increased microbial biomass in grasslands (mean of 8.3%) 21 . As microbial processes play a substantial role in ecosystem models 22,23 , advancing our knowledge of soil microbial responses to moderate, long-term warming is crucial for predicting ecological consequences of future climate changes.
Regions with a Mediterranean climate, characterized by hot, dry summers and cool, wet winters, are particularly vulnerable to global warming [24][25][26] , owing to warming-induced changes such as water scarcity. The Jasper Ridge Global Change Experiment (JRGCE) investigates the impacts of multiple global environmental changes 27,28 , including warming, on a Californian annual grassland experiencing a Mediterranean climate. In this ecosystem, relatively short-term warming (i.e., 2-6 years) had accelerated the flowering and greening of the canopy plants, increased forb production and abundance, stimulated arbuscular mycorrhizal fungi (AMF) hyphal length and total ammonia-oxidizing bacteria (AOB) abundance, and decreased the relative abundance of methane-oxidizing bacteria [27][28][29][30][31][32][33][34][35] . However, it had no impacts on plant biomass, shoot, root, forbs or grass biomass, N cycling process rates, and enzyme activities. At longer time scale (8-9 years), microbial community composition assessed by PLFAs fluctuated greatly over multiple years 36 . Although microbial communities were mostly not responsive to the warming, community composition 1 was changed by warming in 2006. Soil total C concentration also remained unaltered, but the fraction of C derived from microbes, i.e., residues quantified by the microbial amino sugars, significantly decreased in the soils exposed to warming 37 . Those findings are consistent with accumulating evidence that warming often yields different results over the short-and long-terms. For instance, in a temperate deciduous forest, warming had no effect on the composition of the soil bacterial community after 5 or 8 years, but changes were apparent after 20 years 38 . Another 26year warming study conducted in a mid-latitude hardwood forest also showed that soil microbial responses to warming were different over time 39 . Thus, it is essential to gain further insights into long-term ecosystem responses to warming 40 .
Accumulative evidence has noted a decoupling of microbial taxonomy and functional structure 41,42 . Therefore, it is necessary to examine both microbial taxonomy and functional genes. Here, we present data examining the functional and taxonomic composition of soil microbial communities sampled after 14 years of warming at the JRGCE (0.8-1.0°C for 10 years and then 1.5-2.0°C for 4 years). We used Illumina sequencing of 16S rRNA gene amplicons to examine the taxonomic composition of soil bacterial community. We used GeoChip 5.0 to examine functional gene composition derived from bacteria, fungi, and archaea. We also measured plant and soil geochemical variables to examine to what extent changes in taxonomic and functional genes of the soil microbial community in response to warming might be related to changes in plant/soil variables. We aim to test the following hypotheses: (1) Although the warming treatment is moderate, 14 years of warming may be sufficient to induce significant changes in plant communities, microbial communities, and soil geochemical variables; and (2) Microbial responses to warming are different at the taxonomic and functional gene levels.

Plant and soil variables
Warming significantly increased the average soil temperature during the year and month of sampling by 1.1°C at 10-cm depth (Table 1). Perennial forbs biomass (PFB) decreased 49.8% whereas both annual forbs biomass (AFB) and litter biomass (LB) increased 17.1% and 19.9%, respectively, due to warming, though these changes represent a non-significant trend (P > 0.05). Annual grass biomass (AGB), total aboveground biomass (TAGB), total shallow root biomass (TSRB), and fine root biomass (FRB) remained unchanged. Plant standing biomass (PSB), a measure of plant growth consisting of the sum of aboveground biomass (grass and forbs), litter, shallow roots, and fine roots, was not altered by the warming treatment. Soil variables, such as soil moisture, pH, total N, total C, and C:N ratio did not differ between warmed and control plots. Extractable soil NO 3 -N increased 2-fold while extractable soil NH 4 -N decreased 16.7% in response to warming; however, these changes were not significant (P > 0.05).
Taxonomic composition of soil microbial communities Among OTUs obtained from each sample, 7.7% were unclassifiable (not represented in reference sequence databases). Microbial taxonomic composition and alpha-diversity did not differ between the warming treatment plots and the control plots (Supplementary Fig. 1a and Supplementary Table 1), which was supported by results of three nonparametric multivariate statistical tests (P > 0.05, Supplementary Table 2).
The classified OTUs were divided into 2 archaeal phyla and 21 bacterial phyla. The most abundant bacterial phyla were Proteobacteria (27.4%), Acidobacteria (18.7%), Actinobacteria (17.9%), and Verrucomicrobia (7.7%, Supplementary Fig. 2). There were no significant changes in the relative abundances of any phyla or finer taxonomic levels of genus, family, order, or class by warming. A total of 61,908 functional genes, classified into 54,554 bacterial  genes, 4532 fungal genes, 1868 archaeal genes, 400 viral genes,  209 other eukaryotic genes, and 345 unclassified genes, were detected by GeoChip 5.0. Non-metric multidimensional scaling (NMDS) showed that microbial functional gene composition was affected by warming ( Supplementary Fig. 1b). The results of three nonparametric multivariate statistical tests (Adonis, ANOSIM, and MRPP) revealed significant shifts in the functional gene composition of microbial communities in response to warming (P < 0.05, Supplementary Table 2). In contrast, microbial functional diversity and evenness were unchanged (Supplementary Table 1). Functional gene communities of bacteria, archaea, and fungi were changed by warming ( Supplementary Fig. 3), which were also verified by three nonparametric multivariate statistical tests (P < 0.05, Supplementary Table 2).

Microbial functional genes
C cycling. We used the same nonparametric multivariate tests (Adonis, ANOSIM, and MRPP) to determine that the centroids for 19 of 63 functional genes detected by probes related to C cycling were different (unadjusted P < 0.05, Table 2) between 'warmed' and 'control' plots, indicating a significant shift in the composition of C cycling genes in response to warming. These genes include pulA encoding pullulanase, amyX encoding amylopullulanase, nplT encoding neopullulanase, ara encoding arabinofuranosidase, xylA encoding xylosidase, exoglucanase gene, ligninase gene, and endochitinase gene, which were significantly altered by the warming treatment (Table 2). However, these differences were not significant after correcting error rates for multiple comparisons (adjusted P > 0.05, Table 2).
Only seven C cycling genes were significantly changed in abundance ( Fig. 1), as examined by two-tailed paired Student's ttest. Genes exhibiting increased abundances in the warming treatment included amyA encoding alpha-amylase, nplT encoding neopullulanase for starch degradation, and xylA encoding xylose isomerase for hemicellulose degradation. Genes with decreased abundance included mannanase gene for hemicellulose degradation, cellobiase gene for cellulose degradation, acetylglucosaminidase, and chitinase genes for chitin degradation. N cycling. Similar to C cycling, we conducted three nonparametric multivariate statistical tests to examine changes in functional composition of genes related to N cycling. We found that 8 genes (i.e., ureC encoding urease, nifH encoding dinitrogenase reductase, amoA encoding ammonia monooxygenase, nirB encoding nitrite reductase, napA encoding nitrate reductase, nasA encoding nitrate reductase, nirK/nirS encoding nitrite reductase) were altered by warming (unadjusted P < 0.05, Table 2).

DISCUSSION
Our results showed that after 14-year of experimental treatment, warming had no detectable effect on the taxonomic composition and alpha-diversity of soil bacterial and archaeal communities (Supplementary Fig. 1a and Supplementary Table 1). Similarly, a previous study conducted at the same site found that 6 years of warming did not broadly change microbial community composition, as assessed by PLFA markers 36 . The resistance of microbial communities to warming at the taxonomic level might be attributable to several mechanisms. First, microbial communities in Mediterranean-type grassland soils experience considerable yearly and monthly fluctuations in environmental variables, such as temperature and moisture 43 , N supply 44 , and plant growth 36 , which can lead to substantial variations in soil microbial communities 36 . This might explain their weak taxonomic response to warming, as soil microorganisms adapted to frequent environmental changes may better resist further environmental stresses or changes 14,45 . However, treatment effects may play an important role in shaping the microbial communities and soil variables in the long-term 46 , since our data represent the aggregate effects of 14year exposure to continuous warming manipulation despite a single time point sampling. Second, our warming treatment might be not sufficient in magnitude to evoke a significant response in microbial community composition and alpha-diversity because the increase in soil temperature was only 1-2°C, smaller than the magnitude of intra-annual fluctuations in temperature at the time of peak plant biomass 36 . Although warming by infrared heaters typically decreases soil moisture, our study site is quite dry, with average soil moisture in April of 7.9%. Therefore, mild increase in soil temperature of 1-2°C is insufficient to induce notable changes in soil moisture in April (Table 1) 47 , and changes in microbial taxonomic compositions 48,49 . Third, the soil at our study site was warmed 1.0°C for a decade then 1.5-2.0°C for an additional 4 years before sampling. It is possible that 4 years were not long enough to invoke the critical temperature threshold for soil warming that could induce structural changes in taxonomy. Fourth, microbial composition is substantially impacted by plant and soil geochemical variables, such as plant standing biomass, soil organic carbon, pH, soil moisture, and plant coverage [50][51][52] , but these variables did not significantly change with warming in our study. Therefore, our results showed that a 1-2°C elevation in soil temperature, even applied over 14 years, did not have a detectable effect on the overall composition of the soil microbial community. Although fungal taxonomy was not analyzed in this study, we used GeoChip 5.0 to examine functional gene composition derived from bacteria, fungi, and archaea (Supplementary Fig. 3 and Supplementary Table 2). Bacteria usually comprise the majority of microbial biomass in grassland soils 53 , but there is accumulating evidence highlighting the importance of fungal functional roles for grassland ecosystem processes 38,54 . Our results showed that the compositions of archaeal, bacterial, and fungal functional genes were significantly affected by warming. Similarly, in a grassland study of Inner Mongolia, China, functional gene communities of bacteria, archaea, and fungi were changed by light-intensity livestock grazing while bacterial community composition remained unchanged 55 . Bacterial, archaeal, or fungal organisms typically encompass thousands of functional genes, showing a high dimensionality of the functional profile. Most functional genes are not as conserved as the taxonomic biomarker genes (e.g., 16S rRNA gene for bacteria and archaea, and 18S rRNA gene for fungi). Therefore, functional genes are influenced by environmental conditions, leading to high sensitivity to environmental changes.
However, no significant changes were detected in the composition of C cycling genes after correcting error rates for multiple comparisons using Benjamini-Hochberg correction    We found that genes related to labile C degradation pathways increased in warmed plots, whereas those related to recalcitrant C degradation decreased (Fig. 1). This is consistent with results found for other ecosystems, which showed that microbial potential for the metabolization of labile soil C increased in response to warming while that for recalcitrant soil C decreased 62,64 . In addition, a simulated warming study in the alpine grassland of Tibetan plateau showed that relative abundances of catabolic genes associated with more recalcitrant C substrates decreased after warming 60 . A recent study found that 2-year warming increased the total abundance and functional capacities of all potential recalcitrant decomposers in the Alaska tundra 65 . Therefore, the responses of ecosystems to climate warming are complicated, which vary across ecosystem types and the duration of field experiments. Changes in C functional genes might be translated into ecosystems functions, as DNA-based abundances of those genes have been used for assessing CO 2 efflux 66,67 . Warming increased the relative abundance of nifH gene (Fig. 2), suggesting that warming might enhance microbial N 2 -fixation capacity. This result corroborates those from two long-term warming studies conducted in a tall-grass prairie ecosystem in  the Great Plains of Central Oklahoma 63 and a tundra region in Alaska, USA 68 , implying higher N 2 -fixation capacity in response to soil warming. In contrast, warming decreased the relative abundances of many other functional genes related to N cycling. Among those, warming decreased three denitrification genes (narG, nirS, and nosZ), which translated to a 17.7% lower, albeit insignificant, potential for denitrification in our study (Supplementary Table 3). However, potential nitrification rates remained unchanged (Supplementary Table 3). Quantitative linkages between functional gene abundances and process rates have been reported for amoA gene encoding ammonia monooxygenase subunit A [69][70][71] , which mediates the first and rate-limiting step of nitrification in most grasslands. Accordingly, no changes in amoA genes were observed here. In contrast, there was a significant decrease in the relative abundance of hao gene encoding hydroxylamine oxygenase, which mediates the subsequent step of oxidizing ammonia to nitrite through the intermediate hydroxylamine (Fig. 2). There are possible reasons for inconsistency between gene abundances and process measurements. First, we only detected 46 hao genes, far fewer than the 573 amoA genes we report, which suggests a higher degree of functional redundancy for amoA and the process catalyzed by this gene than that for hao. Second, potential nitrification rates may not be indicative of in situ nitrification rates since the environmental conditions are different from those in the field. Third, a previous study conducted at the same site has shown that warming effects on the potential for ammonium oxidation, nitrite oxidation, and denitrification are only detectable when warming has caused changes larger than 30% 33 . In addition, it was shown that long-term warming could increase rates of N cycling 8,62,63,71 .
In this study, we found the reverse situation (i.e., that the relative abundances of most of the key genes involved in N cycling decreased; Fig. 2). This calls for further study, because decreased potentials for key functions involved in N cycling may have important consequences on N cycling in the long term.
Our results indicate that while microbial functional gene abundances and composition are remarkably sensitive to a subtle disturbance of moderate warming, microbial taxonomic composition is remarkably resilient (Supplementary Fig. 1 and Supplementary Table 2). One possible explanation is that functional interchange across different taxa (horizontal gene transfer) and adaptive loss or gain of function through evolution (vertical gene transfer) can collectively result in the conspicuous decoupling of microbial functions from phylogenetic position 72,73 . In addition, we calculated the effective size of warming treatment as the percentages of relative abundance changes: (X t −X c )/X c × 100%, where X t is the mean relative abundance of gene/mean relative abundance of OTU in the treatment, and X c is the mean relative abundance of gene/mean relative abundance of OTU in the control. The smallest detectable effective size for functional gene abundances was 2.5%, meaning that all changes <2.5% were insignificant. In contrast, no significant change was observed for any OTUs, which could be attributed to the finding that 16S rRNA gene amplicon sequencing has lower sensitivity than Geo-Chip 55,64 . However, many questions remain unaddressed. How do environmental conditions differentially interact with taxonomic composition and functional genes? What drives the taxonomic variation within individual functional groups? Can the concept of functional redundancy be used to explain the decoupling of microbial taxonomy and functional genes? In addition, technologies may also play a role. If one performs a metagenomic analysis and extracts taxonomic and functional gene information from the same dataset, will stronger association between microbial taxonomy and functional genes be observed? Future endeavors are needed to tackle those questions.
In summary, here we examine long-term warming effects on soil microbial communities in a Mediterranean-type grassland. We found that warming-induced significant changes in functional gene abundances but had no detectable effect on taxonomic composition of the soil microbial community, suggesting that a functional gene-centric approach to measure and explain responses of microbial communities to environmental changes was more sensitive. Several functional genes involved in the degradation of labile C increased and recalcitrant C decreased by warming, whereas many functional genes involved in N cycling decreased in response to warming. This highlights the importance of taking account of microbial functional potentials in examining ecological consequences of long-term warming. However, an important caveat to bear in mind is that this study lacks much of the insight into interactions that exist between community composition, functional gene expression or activity, and soil processes since there were no measures of soil processes or assessment of gene expression. Future studies should include metatranscriptomic, metaproteomic, or enzymatic activity Fig. 1 The warming effect on carbon (C) degradation genes. The complexity of C is arranged in the order from labile to recalcitrant. Mean abundances of C degradation genes are compared between warmed and control samples, which were calculated as (warming −control)/control × 100%. Error bars represent standard error (n = 4). The differences between warming and control were analyzed by two-tailed paired Student's t-tests (n = 4). **P < 0.05, ***P < 0.01. Fig. 2 The warming effect on nitrogen (N) cycle genes. The percentages in brackets indicate changes in mean abundances of functional genes between warmed and control samples, with red and blue colors representing increase and decrease, respectively. Black color represents no significant changes in the abundances by warming. Gray color represents genes which are not targeted by the version of GeoChip used here, not detected or not applicable. The significant differences between warmed and control sites were analyzed by two-tailed paired Student's t-tests (n = 4). *P < 0.10, **P < 0.05, ***P < 0.01. approaches to generate a better understanding about the warming effects on soil microbial communities.

Site description and experimental design
We conducted this study in an annual grassland ecosystem at the JRGCE site in the San Francisco Bay area, which is located in the eastern foothills of the Santa Cruz Mountains (37°40′ N, 122°22′ W, elevation 150 m), CA, USA. In this Mediterranean-climate ecosystem, plants germinate in October or November and senesce during May and June. The plant community is typical of a Californian annual grassland, and is dominated by annual grasses Avena barbata, Avena fatua, Bromus hordeaceus, and Lolium multiflorum 36,37 , and annual forbs Geranium dissectum, Erodium botrys, and Crepis vesicaria.
This warming experiment is a part of a multi-factor global climate change experiment initiated in 1998, which also involves manipulation of atmospheric CO 2 , N deposition, and precipitation 74 , and an additional fire treatment. The experiment comprised 32 circular grassland plots (2-m diameter plots), each being divided into four quadrants of 0.78 m 2 , resulting in a total of 128 subplots. Every quadrant was separated by a fiberglass barrier (0.5 m belowground) to discourage roots and resources from escaping the quadrants and by a mesh (0.5 m aboveground) to prevent plants, seeds, and litter from crossing quadrant boundaries. Only control subplots (n = 4) and continuously warming subplots (n = 4) were used in the present study; in those subplots, all other global change factors were at ambient levels. Warming was achieved with an array of four overhead infrared heaters (165 × 15 cm; Kalglo Electronics, Bethlehem, PA), which were suspended 1.5 m above the ground in each warmed plot (each heater centered over one quadrant). Temperature was increased by 100 W heaters, resulting in warming of 0.8-1.0°C at the soil surface from 1998 to the end of the 2008-2009 growing season. Then temperature was increased by 250 W heaters, resulting in a warming of 1.5-2.0°C at the soil surface since the beginning of the 2009-2010 growing season. The warming level was within the range of estimates for the second half of this century 75 . Since the spectrum of red: far-red is 622-800 nm, infrared heaters may not have effects on red: far-red received by canopy. Dummy heater hoods were also implemented in control plots to reproduce any shading effects without warming. Soil temperature was measured hourly in each quadrant using thermocouples buried at the depth of 10 cm 33 , and hourly soil temperatures were averaged over the entire month or over the entire year (annual mean).
Soil samples were collected in control and warmed plots on April 26-27, 2012 (i.e., 14 years after the warming treatment started), corresponding to the peak period of plant growth in the four control subplots and in the four warming subplots. One soil core (5 cm in width and 7 cm in depth) was taken from each quadrant. After removing visible plant roots and stones, soil samples were thoroughly mixed and sieved through a 2 mm mesh. Then soil samples were divided into three parts: one part was stored for a few days at 4°C before conducting nitrifying and denitrifying enzyme activity assays, a second part was stored at −20°C prior to a range of soil geochemical measurements, and a third part was stored at −80°C prior to DNA extraction.

Measurements of plant variables
As described previously 29,36 , aboveground plant material was collected on May 5, 2012 from a 141 cm 2 area of each quadrant. Total aboveground biomass (TAGB) was a sum of the biomass of individual plant species. The biomass of individual plant species was combined into functional groups, including annual grass biomass (AGB), annual forbs biomass (AFB), perennial forbs biomass (PFB), and litter biomass (LB). Fine root biomass (FRB) was determined by separating roots out of soil cores taken in the same area of the aboveground biomass harvest. Total shallow root biomass (TSRB) was measured by collecting roots at the depth of 0-15 cm in soil cores. All biomass was oven-dried (70°C) for 72 h before weighing 29,76 .

Measurements of soil geochemical variables and soil microbial activities
Soil moisture of each sample was determined by drying 10 g of fresh soil for a week (or until no change in mass) at 105°C. Soil pH was measured by suspending 5 g soil in 10 ml distilled water and measuring with the glass combination electrode of an Accumet pH meter (Thermo Fisher Scientific Inc., Waltham, Massachusetts, USA). Total soil C (TC) and total soil N (TN) were measured by dry combustion analysis 77  ) under optimal conditions. We incubated sterile Erlenmeyer flasks containing 45 ml of 1 mM phosphate buffer (pH 7.4), 0.04 ml of 0.25 M (NH 4 ) 2 SO 4 , and fresh soil sub-samples (5 g of equivalent dry soil) at 25°C on a shaker with constant agitation at 150 rpm 78 . The 5 ml solution was sampled every 2 h for 10 h and stored at −20°C until analysis of nitrite and nitrate concentrations in a DIONEX ion chromatograph system (Sunnyvale, CA, USA) equipped with an AS11-HC analytical column 79 . Potential nitrification rates were determined from the slope of the linear regression of nitrite plus nitrate concentration on time.
Denitrifying enzyme activity was measured as the nitrous oxide (N 2 O) production rate by soil samples amended with NO 3 − and labile C 80 . Fresh soil (5 g of equivalent dry soil) was placed in a 150 ml plasma flask containing 5 mg glucose, 5 mg glutamic acid, and 0.8 mg KNO 3 . We inhibited N 2 O reductase with a 90:10 N 2 :acetylene mixture. N 2 O concentrations were analyzed every 2 h for 8 h by an Agilent P200 Micro gas chromatograph (Agilent Technologies, Palo Alto, CA, USA) equipped with an electron capture detector 33 . Potential denitrification rates were determined from the slope of the linear regression of N 2 O concentration on time.
DNA extraction, purification, and quantification DNA was extracted from 5 g frozen soil by a freeze-grinding mechanical lysis method as previously described 81 , and then purified using gel electrophoresis with a 0.5% low-melting-point agarose. DNA bands were visualized under UV light and were excised from the gel using a sterile cutter as quickly as possible. Then DNA was eluted in 15 ml double distilled water followed by phenol-chloroform-butanol extraction. DNA was dissolved in sterile, nuclease-free water and stored at −20°C for only a few days before PCR amplification. We assessed DNA purity using the absorbance ratios of A 260 /A 280 (>1.8) and A 260 /A 230 (>1.7). The final DNA concentration was measured with a PicoGreen® method (BMG Labtech, Jena, Germany).

MiSeq sequencing and preprocessing
We used the F515 (5′-GTGCCAGCMGCCGCGG-3′) and R806 (3′-TAATC TWTGGGVHCATCAG-5′) primers to amplify the V4 hypervariable region of bacterial and archaeal 16S rRNA gene. Three replicates of PCR reactions were performed for each sample, which were then pooled together before sequencing on MiSeq platform (Illumina, San Diego, CA) as previously described 82 . Sequencing data were processed and analyzed with the Microbial Ecology Community Pipeline (http://zhoulab5.rccc.ou.edu:8080). Paired raw sequences were identified by unique 12-mer barcodes and merged into longer reads using FLASH 83 . Unqualified sequences were trimmed by Btrim based on quality scores (>20). Sequences were further trimmed to the length of 200-300 bp, and chimera sequences were removed using UCHIME (version 5.2.32) 84 . Reads were assigned to taxonomic units (OTUs) using UPARSE (version usearch 7.0.1001) filtered at 97% nucleotide identity, and a representative sequence of each OTU was generated. To correct for differences in sequencing depth, 28,562 sequences were randomly resampled for each sample before performing the downstream analyses. Sequences observed only once (singletons) were removed to improve data reliability.
Hybridization with GeoChip 5.0 and raw data processing A total of 1 μg DNA was labeled with Cy-3 dye by random priming and hybridized with GeoChip 5.0. The microarray hybridization, scanning, and image processing were conducted as previously described 60,85 . Raw GeoChip data processing was performed with the following steps: (i) remove poor-quality spots with a signal-to-noise ratio (SNR) <2.0, (ii) remove genes detected only in one replicate out of four samples from the same treatment, (iii) perform logarithmic transformation of the signal intensity of each gene, and (iv) normalize the data by dividing the transformed signal intensities of each gene by the mean intensity in the sample.

Statistical analyses
Statistical analyses were performed with R software (v.3.1.0.; R foundation for statistical computing). Functional gene and taxonomic community composition differences were examined by non-metric multidimensional scaling (NMDS) based on the Bray-Curtis distances. Nonparametric multivariate statistical tests (ANOSIM, Adonis, and MRPP) based on different algorithms were performed to examine statistical significance. Specifically, Adonis is an analysis-of-variance analog for multivariate data, which returns an R 2 with a Monte Carlo-permuted P value, ANOSIM is based on rank-order distances and effectively uncovers apparent differences in clusters, while MRPP uses original resemblances and thus is more sensitive to subtle differences across experimental groups. Although they often yield similar results, there are occasions that their results differ. Therefore, it is a popular practice to conduct all three tests simultaneously to show the reliability of the results.
The significance of differences in relative abundances of genes between warmed and control plots were analyzed by two-tailed paired Student's ttests. Normality of the data was verified using Shapiroe-Wilk's tests. Almost all of the genes showed normal distribution and low skewnesses. Only a gene encoding inulinase violates the assumptions of normal distribution, thus the nonparametric analysis, Mann-Whitney U test, was performed. Alpha-diversity indices, including Shannon index (H), Simpsons diversity index (D), and Pielou's evenness (J), were calculated based on the richness and relative abundance data as follows: D ¼ where pi = n i /N, n i is the abundance of the ith OTU (gene), and N is the total abundance of all OTUs (genes) in the sample. S is the species (gene) richness. The P value was adjusted for false discovery rate of 0.05 using the Benjamini-Hochberg method 86 , which was widely used in environmental genomics studies 64,87 .

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.

DATA AVAILABILITY
GeoChip data are available online (www.ncbi.nlm.nih.gov/geo/) with the accession number GSE107168. MiSeq data are available in NCBI SRA database with the accession number SRP126539.