Soil microbial community responses to a decade of warming as revealed by comparative metagenomics

Abstract

Soil microbial communities are extremely complex, being composed of thousands of low-abundance species (<0.1% of total). How such complex communities respond to natural or human-induced fluctuations, including major perturbations such as global climate change, remains poorly understood, severely limiting our predictive ability for soil ecosystem functioning and resilience. In this study, we compared 12 whole-community shotgun metagenomic data sets from a grassland soil in the Midwestern United States, half representing soil that had undergone infrared warming by 2°C for 10 years, which simulated the effects of climate change, and the other half representing the adjacent soil that received no warming and thus, served as controls. Our analyses revealed that the heated communities showed significant shifts in composition and predicted metabolism, and these shifts were community wide as opposed to being attributable to a few taxa. Key metabolic pathways related to carbon turnover, such as cellulose degradation (∼13%) and CO2 production (∼10%), and to nitrogen cycling, including denitrification (∼12%), were enriched under warming, which was consistent with independent physicochemical measurements. These community shifts were interlinked, in part, with higher primary productivity of the aboveground plant communities stimulated by warming, revealing that most of the additional, plant-derived soil carbon was likely respired by microbial activity. Warming also enriched for a higher abundance of sporulation genes and genomes with higher G+C content. Collectively, our results indicate that microbial communities of temperate grassland soils play important roles in mediating feedback responses to climate change and advance the understanding of the molecular mechanisms of community adaptation to environmental perturbations.

FIG 1

Soil community complexity and dominance of sequence-discrete populations. (A) The average coverage, estimated from the portion of nonunique reads (defined as reads with at least one match at the 95% nucleotide identity level; y axis) as a function of the size of subsamples randomly drawn from metagenomes of different habitats (x axis), is shown. The solid lines indicate the fitted model based on subsampling, the empty circles mark the actual size and estimated coverage of the metagenome data sets, and the horizontal dashed line denotes the 95% average coverage level. (B) Eight contig sequences assembled from a control metagenome (C5) were used as references to recruit reads, essentially as described previously (54). The graph shows the identity of each read against the reference sequence (y axes) plotted against the position of the read on the reference sequence (x axes). The histogram on the top represents the read coverage across the length of the contigs; the histogram on the right represents the number of reads recruited per unit of nucleotide identity. Note the genetic discontinuity typically observed in the 95-to-98% nucleotide identity range.

FIG 2

Shifts in taxon abundance and cooccurrence network as effects of warming. (A) Rings represent the average abundances (from six replicate samples) of phyla that made up at least 1% of the whole community; phyla that were significantly different in abundance between heated and control samples are marked by asterisks (P < 0.05, two-tailed paired t test). SPAM is a candidate phylum (55). (B) PCA biplot of the phylum abundance values separated heated from control samples. The contribution of the phyla to each principal component is represented by the arrows. (C) Cooccurrence network, based on Pearson correlation analysis, of the relative abundance of genera in the 12 soil samples (only genera with a correlation coefficient of >0.7 and P value of <0.01 are shown). Each node represents a genus and is color coded for the phylum the genus is assigned to; the size of the node is proportional to the average relative abundance of the genus across the 12 samples. Each line represents a significant correlation between the two genera it connects and is color coded for positive or negative correlation.

FIG 3

Changes in relative abundance of pathways as an effect of warming. The heat map on the left represents changes in the abundance of different pathways (rows) for each pair of samples (columns), color coded based on the magnitude of the change (see scale on the top left). For selected pathways related to the emission of greenhouse gases, the relative abundances of the individual genes that constitute the pathways are shown on the right (small heat maps; rows represent samples, and columns represent genes). In this case, changes in abundance represent deviations from the average abundance of the gene in all 12 samples, are color coded based on the magnitude of the difference (see scale on the bottom right), and are generally consistent with the results for the whole pathway. The results of physicochemical measurements are represented by box plots on the top right. The vertical lines at ratio 1 indicate no change between heated and control samples; the medians of six paired replicate samples are marked by the red bars; the first and third quartiles are represented by the left and right boundaries of the boxes, respectively; the left and right whiskers represent the 1.5 interquartile range; outliers are marked by red asterisks. Abbreviations: SC, soil carbon; MBC, microbial carbon; LPCM1 or -2, labile pool 1 or 2 of carbon, microbial; RspR, respiration rate; SN, soil nitrogen, MBN, microbial nitrogen, LPNM, labile pool of nitrogen, microbial; s- or pMMO, soluble or particulate methane monooxygenase, respectively.

FIG 4

Changes in pathway abundance are community wide and not attributable to only a few taxa. Representative sequences from all OTUs (or clades) of a specific gene (in this case, a CO₂ dehydrogenase, CD_D) were analyzed to produce the distance-based phylogenetic tree shown. The pie charts at the tips of the tree represent the percentage of heated versus control reads that made up each OTU, and the size of each chart is proportional to the number of reads in the OTU; only OTUs with at least 50 reads are shown, for simplicity. Note that no OTU was heat or control specific and about ∼60% of the pie charts had a higher number of heated than of control reads, revealing that many distinct taxa are responsible for the higher abundance of CO₂ dehydrogenase in heated metagenomes. This is also evident in the graph shown at the top (inset). In the circle graph, each circle represents an OTU; the x coordinate represents the percentage of the total reads of the OTU that are control reads, and the y coordinate represents a random value for visualization purposes. Note that more OTUs have less than 50% control reads than have more than 50% control reads. The histogram at the top shows the distribution of the percentages of control reads in all OTUs, with a fitted Gaussian curve shown by the solid line.