|Home | About | Journals | Submit | Contact Us | Français|
The soil ecosystem is critical for human health, affecting aspects of the environment from key agricultural and edaphic parameters to critical influence on climate change. Soil has more unknown biodiversity than any other ecosystem. We have applied diverse DNA extraction methods coupled with high throughput pyrosequencing to explore 4.88 × 109bp of metagenomic sequence data from the longest continually studied soil environment (Park Grass experiment at Rothamsted Research in the UK). Results emphasize important DNA extraction biases and unexpectedly low seasonal and vertical soil metagenomic functional class variations. Clustering-based subsystems and carbohydrate metabolism had the largest quantity of annotated reads assigned although <50% of reads were assigned at an E value cutoff of 10−5. In addition, with the more detailed subsystems, cAMP signaling in bacteria (3.24±0.27% of the annotated reads) and the Ton and Tol transport systems (1.69±0.11%) were relatively highly represented. The most highly represented genome from the database was that for a Bradyrhizobium species. The metagenomic variance created by integrating natural and methodological fluctuations represents a global picture of the Rothamsted soil metagenome that can be used for specific questions and future inter-environmental metagenomic comparisons. However, only 1% of annotated sequences correspond to already sequenced genomes at 96% similarity and E values of <10−5, thus, considerable genomic reconstructions efforts still have to be performed.
Microorganisms first appeared more than 3.5 × 109 years ago (Allwood, 2006, 714–8), ~1.5 billon years after the formation of our planet. Genetic flexibility over a vast expanse of geological time has enabled microorganisms to adapt to virtually every conceivable ecosystem on earth (for example, Huber et al., 1989; Pointing et al., 2009; Larose et al., 2010). Among contemporary ecosystems, soil, which is a product of microbial and macrobial life, exhibits the greatest density and phylogenetic diversity per unit volume (Van Elsas et al., 2006; Roesch et al., 2007), with approximately 109 cells per g, comprising a diversity that is estimated to range from thousands to millions of taxa (Torsvik et al., 2002; Gans et al., 2005).
Soil microbial communities are indispensable for the health of our planet; they drive major geochemical cycles (Falkowski, 2001) and help to support healthy plant growth (Ortíz-Castro et al., 2009). Yet, there is still a considerable lack of understanding of the mechanisms of interaction and metabolism that exist among members of the microbial community and their ecosystem. Existing knowledge, concerning the phylogenetic and functional diversity, community metabolic potential and consequences of evolutionary adaptation, is based largely on partial information gained from studies performed on microorganisms that have been cultivated from soil on a small scale or 16S rRNA gene sequences.
A dependence on studies of cultivable organisms may limit our fundamental understanding of the diversity of interactions in this system. The organisms cultured from soil so far, represent a fraction of the soil biota, for example, those amenable to growth in controlled laboratory conditions (Schloss and Handelsman, 2003; Davis et al., 2011). Attempts to apply metagenomic methodology to soil samples have been hampered by extreme technical challenges, such as extracting an unbiased and representational sample of genetic material from organisms with very different cell membranes and accessible DNA (Handelsman et al., 1998; Ginolhac et al., 2004; Demaneche et al., 2008; Rajendhran and Gunasekaran, 2008; Delmont et al., 2011b, 2011c). This problem is exacerbated by the uneven spatial distribution of microbial communities in soil (Ranjard and Richaume, 2001; Grundmann, 2004). Unlike marine systems, which are generally well mixed and amenable to temporal and biogeographic observations (Gilbert et al. 2009, 2010), soil systems surveys have, despite a wealth of valuable data acquired from hundreds of well-designed experiments and surveys, uncovered only a fraction of the assumed immense microbial diversity of the soil metagenome (for example, Tringe et al., 2005; Roesch et al., 2007; Morales and Holben, 2009). In spite of numerous efforts to study parameters influencing its diversity using cultural-independent approaches (for example, soil pH or nitrogen fertilization, Rousk et al., 2010; Ramirez et al., 2010), data from soil are scarcer than those collected from other commonly encountered ecosystems. The only contemporary published soil metagenome (Tringe et al., 2005) contains just 100 million-bp of DNA, which is potentially a mere millionth of one percent of the genetic material that could be extracted from a g of soil (based on an assumption of 4 million bp per average microbial genome and 109 cells per g of soil). The relative lack of available soil-related sequence data presents an interesting paradox, that the most diverse environment on earth has received the least attention from metagenomic analysis (Vogel et al., 2009) although the first soil metagenome dates from 2005. To redress this balance, we have performed an in-depth investigation of a temperate European ungrazed grassland soil metagenome using pyrosequencing technology.
Building on our previous investigations (Delmont et al., 2011b, 2011c), this study describes an unprecedented effort to characterize the microbial diversity and functional potential of a single soil ecosystem that was found in the Park Grass Experiment at Rothamsted Research; the location of the oldest agricultural experiments in the world run continuously since 1856 (Silvertown et al., 2006). In an attempt to explore this unique environment, almost 5 × 109bp of metagenomic sequence data (Titanium pyrosequencing reads) were produced from soils collected from three depths and at three time points spanning 2 years. To address concerns regarding the influence of DNA extraction technique bias on microbial diversity (Delmont et al., 2011b), we performed 11 different extraction techniques to improve the diversity of the sequenced microbial genomes. The MG-RAST (Meyer et al., 2008) annotated content of the samples were compared with each other, and to samples of two previously reported, non-soil, environments, so as to place the samples in a more global context.
Samples were collected from the untreated control plot (3D) of Park Grass Experiment, Rothamsted Research, Hertfordshire, UK (Silvertown et al., 2006) in February 2009, July 2009 and July 2010. The overall sample handling is outlined in Figure 1. Soil samples from the top 21cm were collected (Delmont et al., 2011b) by sterile manual corers (10cm diameter) in plot 3D at random locations, but not where previous samples had been taken, and were placed in sterile plastic bags, sealed and placed on ice 24h until processing. Previous investigations of this soil demonstrated very little horizontal change in diversity, but measureable changes with depth (Delmont et al., 2011b). Hence, the core samples were fractioned into either seven subsamples as a function of the depth (every 3cm for the direct lysis (described below) and into two depths for the indirect lysis (Delmont et al., 2011c; described below). The aim of this step is to homogenize the quantity of extracted DNA (which decreases with depth) represented in the final pool for each fraction. The different subsamples were then homogenized separately manually by thorough mixing and stored at −20°C for the direct lysis and at 4°C during a maximum period of 1 week for the indirect lysis. To access rhizospheric microbial communities, a soil core (0–21cm) was sieved (0.2mm) and grass roots were extracted. Soil attached to roots was then recovered in a water column. The column helped the physical separation between roots and soil present at its surfaces. The few grams of recovered soil were then mixed before DNA extraction. The metadata for the site and samples are provided in Supplementary Table S1.
Different extraction procedures were used to process the soil samples (Figure 1). We selected DNA extraction methods that use a wide range of approaches to extract and lyse cells. Among the selected methods, some were already known to provide a high DNA yield (for example, BIO1O1), or DNA quality or increased DNA length (in plug lysis), others provided a low yield but could still potentially represent a difficult to access microbial communities. The main goal of this experimental design was to create DNA pools with a large variance in order to uncover a wider range of community members within this soil metagenome at both functional and taxonomical levels.
Utilized one of two bead-beating protocols, (Fast prep MP Bio1O1 Biomedical, Eschwege, Germany) (Griffiths et al., 2000) with 0.5g of soil. This approach was named ‘direct MP Bio101' (F1, J1, and replicates J1a10 and J1b10). In addition, rhizospheric soil from July 2010 was extracted with the same protocol (J1rhizo10) and soil from July 2009 was extracted with another bead-beating method, the MoBio PowerSoil DNA Isolation Kit (Carlsbad, CA, USA) (J7). Several different indirect DNA extraction methods were used by first extracting cells on a Nycodenz gradient gel (density of 1.3) (Bertrand et al., 2005) and then applying one of the following four lyses with the extracted cells: (1) the same bead-beating protocol, called ‘indirect MP Bio101' (replicates F2a and F2b from February 2009); (2) the Nucleospin Tissue kit, named ‘indirect DNA Tissue' (F4 and J4, February and July 2009, respectively); (3) the Gram-positive kit, named ‘indirect Gram positive' (F5); and finally (4) a lysis using agarose plugs called ‘indirect lysis in plug'. (F3 and F6 from 0 to 10cm and 11 to 21cm depths, February 2009, respectively—see Figure 1). Plugs were first transferred in 3ml of G− lysis buffer (1% lauroyl sarcosine, 500mM of EDTA Na2, pH 9.5) with 0.5mgml–1 of lysozyme and incubated at 37°C for 12h. The agarose plugs were then incubated in 3ml of G− lysis buffer with 500μgml–1 of proteinase K at 56°C for 12h, and finally equilibrated in a 10mM Tris (pH 8.0), 1mM EDTA storage buffer). This enzymatic lysis was performed in a stable environment (the agarose plug) and was performed without any physical perturbations (for example, mixing in the tubes resulting in DNA breakage). This method is generally used to provide high quality and long DNA sequences for the construction of fosmid libraries or for genome size. General information about the different DNA extraction yields used is presented in the Table 1.
A minimum of 10μg of DNA were used for each Roche/454 pyrosequencing run on a 454 pyrosequencer (GS FLX Titanium Series Reagents; Roche 454; Shirley, NY, USA). Processing of samples (before sequencing) did not involve prior amplification step. For the direct lysis, equal quantities of DNA extracted from the seven fractions from 0 to 21cm were pooled together. J1a10 and J1b10 correspond to distinct extractions from the same soil core. For the indirect approach corresponding to soil from 0 to 21cm, equal quantities of DNA extracted from the two fractions (0–10 and 11–21cm) were pooled together. For the indirect lysis using the bead-beating protocol (0–21cm, February 2009), two pyrosequencing runs (F2a and F2b) were performed from the same DNA pool (>20 micrograms). The sequence data are publically available (http://www.genomenviron.org/Projects/METASOIL.html).
Artificial duplicates were deleted using cd-hit-454 with default parameters (Niu et al. 2010). Sequences were then annotated on the MG RAST (v.02) online software (Meyer et al., 2008). Reads were distributed into different metabolic subsystems. Similarity search between pyrosequencing reads and the SEED database (Overbeek et al., 2005) have been processed with a maximum E value of 10−5. All compared distributions were normalized as a function of the number of annotated sequences for each metagenome. Data corresponding to both functional and taxonomical distributions were then statistically analyzed within the STAMP software (Parks and Beiko, 2010). Fisher's exact tests were performed and annotated functions and taxa with P-values <0.05 were considered to be significantly different between the different experiments.
Tests on assembly productivity were performed using Newbler (Margulies et al., 2005). Newbler was run directly from the ‘.sff' files produced by the pyrosequencer using the following parameters: expected depth, 0 (that is, undefined); minimum read length, 20; seed step, 12; seed length, 16; seed count, 1; minimum overlap length, 40; minimum overlap identity, 90% alignment identity score, 2; alignment difference score, −3. The minimum read length in the data set was 40. Each deeply sequenced data set was assembled separately using the 454 GS de novo assembler software (Newbler v2.0.00.22), and all contigs were used for subsequent analysis. In addition, MetaGeneMark (version 2.7d using the parameter file for metagenome gene prediction version 1; http://exon.gatech.edu/metagenome/Prediction/) have been used to search genes from the 100 largest contigs, and the 1006 genes predicted were analyzed via MG-RAST.
A total of 13 pyrosequencing runs were performed with DNA extracted from the Rothamsted Research (Park Grass) site. Grassland soil samples were taken at different depths and three different time points over 1.5 years. DNA was extracted from the samples using six different DNA extraction protocols (see Materials and methods and Figure 1). Two samples were sequenced in duplicate (J1a10 and J1b10, and F2a and F2b) to explore the reproducibility of the metagenomic profile. A total of 12575129 reads were generated (length average of 385.9±31.8bp) and 34.5 (±3.3)% of them were annotated with the MG-RAST online server (E value <10−5) (Meyer et al., 2008). On the basis of the protein database used by MG-RAST, 88.64 (±1.44)% of these annotated sequences had closest homology to a protein found in Bacteria, 0.91 (±0.23)% to Eukarya, and 1.41 (±0.16)% to Archaea. Thus, almost 9% of annotated sequences were not classified at the domain level. All the annotated reads were compared with SEED-NR, FIGFams for functional assignments and then used in subsystem reconstructions. The closest matched gene was the source of information about the functional (metabolic) subsystem that the read was binned into and about the ‘taxa' represented by this read. Therefore, the taxa cited here correspond to the genomes in the database that best matched the given read as long as the E value was <10−5. Major functions and taxa identified can be found in Supplementary Tables S2 through S4.
Functional differences between the 13 datasets generated from Rothamsted, two datasets from other soils, and one from an aquatic environment were derived by exploring the relative number of reads associated with the 835 functional subsystems detected at least in one metagenome (Figure 2). Bootstrap values are provided. The method of DNA extraction correlates with sample grouping. Samples, F1, J1, J1.a and J1.b were directly extracted using the MP Bio101 kit. Sample J7, which lies in the same general group was extracted directly with MoBio Powersoil kit. The sample from the application of direct MP Bio101 on the rhizosphere soil (J2) is closely associated with this group. The bootstrap values are not particularly high within this group. Three sample pairs on the other hand had significant bootstrap values (>90%) grouping them apart from the other samples: (1) the replicate samples from the application of MP Bio101 to the cells first removed from soil via the Nycodenz gradient (F2.a and F2.b); (2) the two depth samples extracted by indirect lysis in agarose plugs (F3 and F6); and (3) the two samples from different seasons extracted after Nycodenz by use of the DNA tissue kit (F4 and J4). In order to assess the statistical likelihood of the subsystem distribution differences between samples, STAMP software (Parks and Beiko, 2010) based on a bootstrap approach using Fisher's exact tests was applied to the MG-RAST (subsystem functional level 3) outputs. This approach determined what percentage of the 835 subsystems were significantly (at 95% CI) different between any pair-wise comparison. Replicate runs (F2a/F2b and J1a10/J1b10) had between 7.3% and 7.7% dissimilar subsystems, and seasonal variations had 8.6% and 11.7% dissimilar subsystems for direct (F1/J1) and indirect (F4/J4) extractions, respectively. When different lysis methods (for example, F4, F5 and F6) were applied to the bacterial cells removed by Nycodenz gradient gel before DNA extraction, significant differences in subsystem distributions (16.9–39.8% dissimilar subsystems) were observed at the 95% CI. Using sequences corresponding to communities extracted from two distinct horizons (0–10cm: F3 and 11 to 20cm: F6), 27.01% of the detected functional subsystems possessed statistically different distributions. Two types of geographical comparisons were made. One was between Rothamsted soil and soil from Italy (Vallombrosa forest soil, defined as a Cambric Umbrisol) extracted with the same method (MP Bio101) in our laboratory (the related Italian soil metagenome is represented by approximately 100000 sequences) and these two soils had 14.1% dissimilar subsystems. The second was Rothamsted sequences compared with those from Puerto Rico (located in the Luquillo experimental forest and defined as a tropical rain forest soil, Metagenome ID of 4446153.3 on MG RAST, one million reads), which were extracted and sequenced elsewhere. They had between 30.98% and 33.13% dissimilar subsystems. The most extreme comparison was between Rothamsted soil and the Sargasso Sea (72% dissimilar subsystems) as indicated also by the distance in Figure 2.
Among the major (29) metabolic classes, clustering-based subsystems (CBSS) and carbohydrate metabolism had the largest quantity of annotated reads assigned (Figure 3). Virulence and amino acid and derivatives were next in prevalence (Figure 3). The cluster-based subsystems contain such functions as proteosomes, ribosomes and recombination-related clusters. The virulence subsystem contains diverse functions also, such as resistance to antibiotics and toxic compounds, and pathogenicity islands. Some subsystems were relatively minor such as photosynthesis, prophage, dormancy and sporulation (Figure 3). Although there was significant (at the 95% CI) differences in the distribution of reads in some (from about 7 to 40%) of the different metabolic subsystems from the different pyrosequencing runs of DNA extracted from Rothamsted soil, the s.d. around the mean of all of the pyrosequencing runs varied between 2% and 50%, with the higher variance for the metabolic classes with relative few assigned reads (for example, macromolecular synthesis; error bars in Figure 3). When comparing the assigned reads at a finer functional subsystem classification within MG-RAST, the most prevalent subsystem (out of the 835 different categories) in the soil was the cAMP signaling in bacteria with 3.24%±0.27% of the annotated reads (Supplementary Table S2). The next most prevalent subsystem was the Ton and Tol transport systems at 1.69%±0.11% of the annotated reads (Supplementary Table S2). These prevalent systems varied <10% between DNA extraction pools except for the distribution of CO2 uptake carboxysome-related genes, which varied from 0.56% in F3 to 1.43% in F4, which represents an increase of 60.7% (average for the 13 pyrosequencing runs was 0.70%±0.42% of reads).
The 56%±4.4% of annotated protein sequences showed closest homology to a total of 1214 unique taxa using the taxonomic annotation of functional SEED subsystems. The most dominant putative taxon was Solibacter usitatus (6.72%±0.29% of annotated reads). Other taxa with relatively high number of assigned reads were Blastopirellula marina (4.96%±2.88%), Bradyrhizobium japonicum (4.89%±0.635%) and Acidobacteria bacterium (3.64%±0.94%) (Supplementary Table S3). The legitimacy of the read assignment at an E value of 10−5 cutoff is provided in part by the E value distribution of the different reads assigned to the reference genome. In the case of B. japonicum, the majority of assigned reads had E values <10−30 and in the case of B. marina, the E values were in general >10−30. With the taxonomic classification of functional gene fragments, it is possible to use all annotated reads to determine community structure (Figure 4); however, it is also possible to use 16S rRNA sequences to determine the community structure, although the number of reads is considerably less than for the SEED annotation. Three different databases accessible within the MG-RAST platform and used with the MG-RAST software were used to determine community structure with s.d. calculated from the variance of the 13 different pyrosequencing runs (Figure 4). Although there is a general agreement: alpha-, beta- and gammaproteobacteria and Actinobacteria dominate all four methods, there are some important differences in the relative number of reads in different classifications. For example, the Silva SSU database 94 has a much higher percentage of reads in Flavobacteria than the other systems (Figure 4). In order to use functional genes other than 16S rRNA for taxa or at least genera identification, a more accurate and limited analysis constrained the similarity at 96% or better (still with an E value of 10−5). When this was performed using SEED, only 0.35%±0.09% of the total reads (or about 1% of the annotated reads) were used to identify bacterial taxa (Supplementary Table S4). The most abundant taxa identified from Rothamsted soil were members of the Bradyrhizobium, Rhodopseudomonas and Nitrobacter genera (Alphaproteobacteria); the Solibacter and Acidobacteria genera (Acidobacteria) and Pseudomonas (Gammaproteobacteria) and Burkholderia (Betaproteobacteria) genera. B. marina was no longer associated with any of the reads.
Sequence data were assembled to provide a metric describing the depth of sequencing applied to the community metagenome; in part, this was used to estimate the minimum quantity of sequencing required to completely sequence all the members of the soil microbial community. The extreme minimum could be considered as the quantity of sequences where no singleton is left unassembled, even if practically this minimum is insufficient to assemble all the genomes. A total of 10 random read subsamples of increasing metagenome size (read quantity) were run through the Newbler assembler (Supplementary Table S5). No attempt was made here to optimize the assembly process. The 10 subsamples ranged in size from 1257242 to 12572342 reads (with 487554794 to 4874169257 total number of bases) and produced from 7478 to 266600 contigs (Supplementary Table S5). The largest contig size increased from 6361 to 22645bp with three times as many reads but then decreased and leveled off at about 15400bp with increasing number of reads (Supplementary Table S5). The fraction of reads that were not included in any contig (‘singletons') fell from roughly 0.93 to 0.76 when increasing the number of reads 10-fold (Figure 5a insert). This data was fitted and extrapolated to the point where no read would be orphaned. This extrapolation was at about 400 million 454 reads (average of 386bp in length) with the 95% confidence interval stretching from <200 million reads to almost 1400 million reads (Figure 5a). The maximum contig length did not continue to increase with increasing read number, but the number of reads per contig did develop two general trends (Figure 5b). The denser trend has a slope represented by a contig coverage of about 30 × (when the assembler needs/uses 30 × to build the contigs) and the smaller trend has a contig coverage of about 4.5 × . Contigs from these two trends were selected, broken in coding sequences by MetaGeneMark and then annotated using MG RAST. Globally, the trend corresponding to coverage of 30 × possessed more sequences related to Firmicutes (10.99%) and Verrucomicrobia (21.85%). In contrast, the trend corresponding to low coverage assembled contigs (4.5 × coverage) possessed a majority of sequences related to Proteobacteria (66.06%). Independent of the two observed trends, the 100 largest contigs created from the entire sequence pool were also annotated by MG RAST and in general the relative proportion of different functional and phylogenic classes (stars in Figures 3 and and4)4) were similar to that for the sequences directly with some exceptions. There were fewer virulence subsystem hits and significantly more fatty acids and protein metabolism hits.
Soil is one of the most diverse environments on earth and the depth of the microbial diversity is still poorly understood. High throughput sequencing technologies, coupled with appropriate DNA extraction methods, provide a means to explore the soil ecosystem with an unprecedented level of detail (Vogel et al., 2009). In this study, pyrosequencing from 13 samples generated nearly 5 × 109bp of sequence data with average read size of 386bp. Three key parameters were varied: soil depth, sample collection season and DNA extraction method. Sequence samples were annotated with the MG-RAST online server, revealing broad functional (835 of 878 possible functional subsystems) and taxonomic (detection of 1214 putative taxa) diversity in the Rothamsted Park Grass soil metagenome.
The most abundant functional subsystems in the Rothamsted soil seemed to be related to microbial cAMP signaling and Ton and Tol transport (Supplementary Table S2). The same subsystems were prevalent in metagenomes in soil at Waseca farm, in Puerto Rico and Italy. These trends in soil functional content are robust enough to be observed on a global scale. cAMP is an important secondary messenger in Eukarya and Bacteria. cAMP is a universal cell energy/metabolism regulator as well as being involved with cell–cell signaling. Soil bacteria might have to deal with frequently fluctuating substrate levels so that they would need extra regulation rather than interacting with plants. Interestingly, since cAMP is also a subversion mechanism, some bacterial pathogens might also subvert plant cAMP production for their own benefit, through injection of adenylate cyclase and/or various toxins that alter adenylate cyclase levels (adenylate cyclase is essential to the production of cAMP) (Akhter et al., 2008; Agarwal and Bishai, 2009). Iron is an essential element for most organisms (Weinberg, 1984), but can be a limiting reagent for life (often in oceans, Boyd et al., 2007) owing to its insolubility in aerobic environments at neutral pH. In response to this stress, some bacteria possess high-affinity transport systems (Crosa et al., 1997) and generate high-affinity siderophores that complex extracellular iron to optimize its acquisition. The presence of Ton-related proteins in the soil is likely due to TonB, an energy-dependent cell envelope protein that assists iron uptake through accommodation of ferric siderophores, too large to cross porins, through the outer membrane (Klebba et al., 2003).
MG-RAST annotation also revealed the presence of several highly abundant CBSS. These are groups of functionally coupled genes (genes found proximal to each other in the genomes of diverse taxa) whose functional attributes are not well understood. The relatively high abundance of these subsystems across all Park Grass samples, as well as the other sequenced soils, suggests that they have key roles in soil ecosystems across the globe, and should be explored in future research efforts to understand the composition of soil ecosystems. The CBSS-258594.1.peg.3339, CBSS-269799.3.peg.2220, CBSS-83332.1.peg.3803, CBSS-249196.1.peg.364 (Supplementary Table S2) are thought to be a galactoglycan biosynthesis, a molybdenum oxidoreductase, a PKS-related, and a fatty acid metabolism subsystem, respectively.
The comparison of the runs corresponding to the same DNA sample (F2a/F2b) provided important information about the reproducibility of pyrosequence generation in highly biodiverse environments. The Fisher's exact test operated by the STAMP software did identify some functions (about 7%) and taxa that varied significantly (at the 95% CI) between replicates. The lower P-value was on the order of 10−7 when comparing F2a and F2b at the functional level, so some comparisons between seasons and depths were possible. On the basis of these observations, functional comparisons having at most a minimum P-value of 10−8 (cutoff based on the observed technological reproducibility) were considered to have distributions that varied significantly. Unfortunately, the technological reproducibility is not the only limit for robust metagenomic comparisons. Even if a stringent P-value is used, the DNA extraction approach influenced the experimental conclusions. When comparing the seasonal effect by using two different extraction approaches (direct:F1/J1 and indirect F4/J4), some differences in relative predominance of different subsystems were found. On the basis of the comparison of F1 and J1, sequences related to the type 4 secretion and conjugative transfer and cellulosome subsystems are more represented in February (P-value of 10−8 in the two cases). When comparing F4 and J4, the cellulosome subsystem is still detected more in February (P-value<10−15) but the type 4 secretion and conjugative transfer is not. In contrast, sequences related to bacterial cAMP signaling are more present in July (P-value of 10−12), but only when comparing F4 and J4. Thus, only sequences related to cellulosome dominated one season's metagenome independent of the extraction method applied. Major environmental difference between the two studied seasons was temperature (from 6°C in February to 16.6°C in July). In addition, snow lay on the ground for weeks in February of the same year, thus limiting active grass growth. As a consequence, soluble root exudates were possibly in short supply during this relatively cold period and cellulosome from root residues would be the main source of carbon and energy supporting soil microbial communities.
On the other hand, depth had more effect with sequences related to genes involved in bacterial chemotaxis, Ton and Tol transport systems, flagellum mechanism, D-ribose and L-Arabinose utilization represented more in the surface sample (0–10cm) and sequences related to selenocysteine metabolism and tRNA aminoacylation represented more at depth (11–20cm). However, these results were generated using only one DNA extraction method. In comparison to depth and seasonal variables, the extraction method was able to influence functional distributions (Figure 2), especially when using methods with striking differences in cell lysis (for example, Gram positive kit versus in agarose plug lysis or DNA tissue). Thus, the stringency of lysis appears to be a crucial step for soil metagenomic analysis, confirming previous results with RISA and phylogenetic microarray analyses (Delmont et al., 2011b).
In addition, when studying the distribution of sequences based on their G+C%, clear variations were found among the different runs. Direct lysis versus indirect lysis had more impact on the G+C% profile than any other variable. The indirect lysis provided more sequences possessing a higher G+C ratio (from 60% to 72%), whereas the direct lysis had a more even distribution with more sequences in the 50 to 58 G+C% range (Supplementary Figure S1). Both metagenomic s.d. and G+C% ratio profile fluctuations are limited by the experiments and variables used. However, this effort provides both significant soil metagenomic sequences and data useful to appreciate methodological differences in microbial community diversity accessibility.
Given the relatively low functional subsystem variations between different soils (Figure 2), soil microbial community metagenomes from Rothamsted, Puerto Rico, Italy and the Waseca farm soil (Tringe et al., 2005) could be compared with metagenomes from oceans and human feces. This comparison might help identify some of the soil ecosystem unique functional attributes. In order to make the comparison, principal component analysis was generated based on the distribution of general functional subsystem classes with metagenomes publically available from these ecosystems (Figure 6). Some general functional classifications appear to be relatively more represented in one ecosystem in comparison with the others. Sequences related to RNA and protein metabolism, photosynthesis, fatty acids and lipids, and macromolecular synthesis are more highly represented in ocean metagenomes. In contrast, phosphorus metabolism and virulence are less represented in ocean metagenomes than in those sequenced for soil and human microbiomes. Sulfur and potassium metabolism, membrane transport, stress response and regulation, and cell signaling are more represented, and nucleosides and nucleotides, and RNA and protein metabolism are less represented in soil metagenomes. In human microbiomes, cell division and cell cycle, DNA and phosphorus metabolism, cell wall and capsule, dormancy and sporulation, and carbohydrates are more represented than in those of oceans and soils (Figure 6). When comparing the taxonomical structure of these metagenomes, Cyanobacteria and Bacteroidetes appear to be more represented in the oceans. In addition, eukaryotic sequences were also detected and represent additional specificities of these metagenomes (Supplementary Figure S2). Actinobacteria, Chloroflexi, Fibrobacteres and Acidobacteria group, Planctomycetes, and Synergistetes are more numerous in soils. Chlorobi, Firmicutes, Spirochaetes, Fusobacteria and the Bacteroidetes Chlorobi group are clearly relatively dominant in human digestive tracts. In contrast, more Proteobacteria are present in oceans and soils. The metagenomes are clearly grouped as a function of the environment based on both general functional and taxonomical distributions. So in spite of important DNA extraction biases and sequencing technology differences (Illumina, Pyrosequencing and Sanger), global metagenomic comparisons are possible and provide unique information about the functional and taxonomical differences of each environment (Delmont et al., 2011a). As an example, sequences related to metabolism of aromatic compounds are more abundant in soils possibly due to the presence of these compounds in this environment. However, additional comparisons, such as qPCR and metatranscriptomics, need to be performed to confirm which taxa and functions are unusually active in soil to gain a better understanding of soil microbial community function.
The relative percentage of orphan reads decreased continually when accumulating pyrosequences. Therefore, an estimate of the number of reads needed to avoid having orphan reads would possibly provide the absolute minimum number of reads needed to sequence the entire soil metagenome. Rarefaction analysis of this sequencing effort (Figure 5) indicated that the equivalent of about 450 Titanium runs would be required to create contigs from all of the soil pyrosequence reads generated. Of course, chimeras might be generated due to the complexity of communities, and a much larger effort would be needed to assemble the soil metagenome, but as new efficient high-throughput sequencing technologies and valuable assembling tools are developed, this goal will become less utopic. Genomes from Proteobacteria might be assembled more rapidly than those from Firmicute or Verrucomicrobia phyla. The presence of regions that limit assembly (for example, insertion sequences regions) and the complexity of diversity among taxa might explain in part the efficiency differences observed between these phyla (4.5 × and 30 × ), but additional experiments are needed to understand the two trends observed in the Figure 5B.
In this study, >12 million reads were generated from the soil of the Rothamsted Research Park Grass experiment. These sequences were generated in 13 separate sequencing runs producing over 4 × 109bp. The results demonstrated both some DNA extraction biases and relatively low seasonal (when comparing February and July months) and vertical soil metagenomic functional class fluctuations. In addition, this approach provided a statistical view of functional distributions in this soil. This metagenomic study increased our knowledge about soil microbial communities at a metagenomic level by integrating both natural and methodological fluctuations. The metagenomic variance so generated represents a global picture of the Rothamsted soil metagenome that can be used for specific questions and future inter-environmental metagenomic comparisons. However, only 34.5% of the reads were assigned to functions and <1% of annotated sequences correspond to already sequenced genomes (at 96% similarity), therefore, many soil microorganisms remain elusive and genome constructions are needed.
TOD was supported by the Rhône-Alpes Region. We would like to thank the French National Research Agency (ANR) for financing Metasoil (Projet ANR-08-GENM-025) and the European Union (7th Framework KBBE-2007-3-3-05) funding for Metaexplore (22625) project.
Supplementary Information accompanies the paper on The ISME Journal website (http://www.nature.com/ismej)