Methane-Oxidizing Bacteria in a California Upland Grassland Soil: Diversity and Response to Simulated Global Change

The impact of individual and multiple, simultaneous global changes on methanotroph community composition was investigated using the Jasper Ridge Global Change Experiment (JRGCE). The JRGCE is located on the Jasper Ridge Biological Preserve, which lies in the eastern foothills of the Santa Cruz Mountains in northern California. The climate, vegetation, and soil parameters, as well as the experimental design, have been described in detail previously (43, 48). In brief, the JRGCE was established in a grassland ecosystem dominated by annual grasses (Avena barbata and Bromus hordeaceus) and forbs (Geranium dissectum and Erodium botrys), growing on a sandstone-derived soil with an average pH of 6.31 ± 0.3. Four global change factors, CO₂ (ambient and 680 ppm), temperature (ambient and ambient plus 80 W m⁻² of thermal radiation), precipitation (ambient and 50% above ambient), and nitrogen deposition (ambient and ambient plus 7 g N m⁻² in the form of calcium nitrate), were applied to different plots in a full factorial design (leading to a total of 16 different treatments). Each treatment was replicated eight times. The treatments were applied as a split-plot design with 32 circular plots, each divided into four 0.78 m² quadrants, separated by solid belowground and mesh aboveground partitions. Infrared heat lamps were suspended over the centers of the warming plots, heating the plants in all quadrants of a plot by 0.8 to 1°C. Atmospheric CO₂ concentrations were elevated with a ring of free-air emitters surrounding the plots. Ambient precipitation events were augmented with drip irrigation and overhead sprinklers; the precipitation treatment increased the average soil moisture from 19.8% to 26.6% (measured at the time of soil sampling). Warming and CO₂ treatments were applied on the whole-plot level, and precipitation and nitrogen treatments were applied on the subplot level. Manipulations started in the autumn of 1998, at the beginning of coastal California's rainy season.

The analysis of microbial communities was initiated in May 2000. Soil cores from all replicate treatments were taken from a depth of 15 cm with a 2.2-cm-diameter corer. Each core was placed in a plastic bag, cooled on ice in the field, and homogenized thoroughly by hand in the laboratory prior to storing at −80°C.

Extraction of DNA from 0.5 g of soil was performed using the Ultra soil DNA extraction kit (MoBio Laboratories, Solana Beach, CA) according to the manufacturer's instructions, with the exception that the final purification step was repeated to increase the purity of the DNA. The DNA was resuspended in a final volume of 50 μl and stored at −80°C. DNA quantification was performed with the PicoGreen assay (Molecular Probes, Eugene, OR) according to the manufacturer's directions. The DNA yield was approximately 5 to 20 ng/μl.

To characterize the methanotrophic diversity, we tested five different primer combinations for their suitability to amplify pmoA gene types in the Jasper Ridge grassland soils. For this preliminary test, we chose soil samples from two plots with elevated CO₂, temperature, precipitation, and nitrogen (plots ID5 and ID60). For each soil and primer combination, one clone library was generated, and we sequenced 15 clones per library. The primer combinations tested were (i) A189F-682R (24), (ii) A189F-650R (10), (iii) A189F-mb661R (12), (iv) A189F-682R (seminested; 650R), and (v) A189F-682R (seminested; mb661R). All clones sequenced from clone libraries generated by use of the A189F-682R primer system were amoA sequence types closely related to the ammonia-oxidizer Nitrosospira multiformis. Clone libraries that were generated based on the A189F-650R and A189F-mb661R primer systems contained some pmoA sequences. However, up to 50% of the randomly selected clones contained nonspecific inserts. In contrast, all clones sequenced from clone libraries generated using the two seminested PCR approaches, A189F-682R (seminested; 650R) and A189F-682R (seminested; mb661R), were pmoA sequence types. Therefore the seminested PCR approach was subsequently used for the study of methanotrophic diversity and for generating community profiles by T-RFLP analysis.

As described above, the amplification of pmoA genes was performed via a seminested PCR approach using the 5′ primer A189 and the 3′ primer A682 (24). The temperature profile (Table 1) was identical to the previously described “touch-down” PCR protocol (28). Aliquots of the first round of PCR (0.25 μl) were used as the template in the second round of PCR using the 5′ primer A189 and the two 3′ primers mb661R and 650R in a multiplex PCR setting (i.e., both reverse primers were present in the same reaction). This approach allowed simultaneous amplification of a broad range of pmoA targets. The reverse primer mb661R was designed for the detection of type I and type II methanotrophs (12), while the reverse primer 650R was designed for the specific detection of putative atmospheric methane oxidizers from the “RA 14” clade (10). Each reaction mixture contained 12.5 μl of MasterAmp PCR premix F (Epicentre Technologies, Madison, WI), 0.5 μM of (each) primer (QIAGEN, Alameda, CA), 1.25 U of Taq DNA polymerase Low DNA (AmpliTaq, Applied Biosystems, Foster City, CA), and 0.25 μl of template DNA. Amplification was performed in a total volume of 25 μl in 0.2-ml reaction tubes, using a DNA Engine thermal cycler (MJ Research, San Francisco, CA). The PCR amplifications of environmental DNA resulted in amplicons of the expected size (approximately 500 bp). The first round reaction and the second round reaction were each performed in triplicate. Aliquots from the first round (three independent reactions in three different tubes) were pooled before going into the second round (which was itself done in triplicate). These final reactions were pooled prior to digestion. Aliquots of the amplicons (5 μl) were checked by electrophoresis on a 1% agarose gel.

PCR products were cloned using a TOPO TA cloning kit (Invitrogen Corp., San Diego, CA) following the protocol of the manufacturer. The preparation of plasmid DNA of randomly selected clones, PCR amplification of cloned inserts, and nonradioactive sequencing were carried out as described previously (28).

The identities of the pmoA gene sequences were confirmed by searching the international sequence databases using the BLAST programs (http://www.ncbi.nlm.nih.gov/BLAST/ ). The currently available database of pmoA gene sequences was integrated within the ARB program package (33), and DNA sequences were analyzed and edited using the alignment tools implemented in ARB. We constructed phylogenetic trees using the maximum likelihood approach (with the default settings), the Fitch-Margoliash approach (using global rearrangement and randomized input order with three jumbles), and the neighbor-joining approach (with the Felsenstein correction) in ARB. The robustness of the tree topology was verified through calculating bootstrap values for the neighbor-joining tree and through comparison of the topology of the trees constructed using the different approaches.

The molecular evolution of the novel pmoA lineages was investigated using the codeml executable of the Phylogenetic Analysis by Maximum Likelihood (PAML) program (58). The input nucleotide files contained a 453-nucleotide portion of all pmoA sequences shown in Fig. 1 (with the exception of sequence “E5FB-b” [AJ579668 ] from the “WB5FH-A” clade, which was too short to include, and “LOPA 12.6” [AF358043 ], “IY-6.48” [AY236518 ], and “RA 14” [AF148521 ], which were added to Fig. 1 during final revisions of the manuscript; in addition, “VIp9” [AY37258 ] was removed from Fig. 1 during final revisions, but was present in the PAML analysis). To reduce the level of sequence divergence to within recommended levels (Z. Yang, personal communication), the PAML analyses were run on the two halves of the pmoA phylogenetic tree separately. One half contained the type I “WB5FH-A,” JR2, and JR3 pmoA clades along with the two Nitrosococcus amoA sequences. The other half contained the type II “RA 14” and JR1 clades as well as Methylocapsa acidiphila. Due to the divergence of the two amoA sequences from the pmoA sequences in members of the class Gammaproteobacteria, the type I side of the tree was still at the limits of acceptable divergence for the PAML program, and so analyses for that side of the tree were also run without the amoA sequences. All sequences were in frame and aligned (using MacClade 4.03 PPC; Sinauer Associates, Inc., Sunderland, MA), and the few ambiguous sites were assigned the nucleotide of their nearest phylogenetic neighbors. The input tree files were created using PAUP 4.0b10 (49), using analysis by distance, neighbor joining with Jukes-Cantor correction, and ties broken randomly. Their topology matched that of the tree (Fig. 1) presented in this paper.

Branch lengths were estimated by the PAML program using the one-ratio model, and then those branch lengths were used as the initial values for branch length estimation in further models performed. In the codeml control file, the majority of parameters were left in their default specifications, with the following exceptions: runmode = 0, seqtype = 1, CodonFreq = 2, Model = 0 or 2, and, for the multiratio models, fix_blength = 1.

The one-ratio model was run to provide an estimation of a single nonsynonymous-to-synonymous substitution (“dN/dS”) ratio for each half of the tree. A series of two-ratio models were then run, to allow the dN/dS ratio of the three novel lineages to vary in turn. Lastly, the dN/dS of each major branch and clade (as denoted in Fig. 2) was allowed to vary simultaneously under the freely varying model, generating maximum likelihood estimates for all dN/dS values across the tree (57).

To test the robustness of the parameter estimates, all analyses were also run on various subsets of the taxa, with little variation in the results; this is consistent with other studies that have shown that codeml is robust to sampling (57, 59). In addition, all analyses were run at least twice to ensure that parameter estimates were likely global rather than local optima.

The likelihood ratio test was used to assess the goodness of fit of the two-ratio models to the data and to compare it with that of the one-ratio model. This allowed us to test whether the dN/dS ratios on the branches leading to the three novel clades were significantly different from the background dN/dS ratio in the remainder of the tree (57).

The creation of terminal restriction fragments (T-RFs) from pmoA genes was carried out as previously described (28). After purification with QIAquick spin columns (QIAGEN, Alameda, CA), approximately 100 ng of the amplicons was digested separately with 20 U of the restriction endonuclease MspI (New England BioLabs, Beverly, MA). The digestions were carried out in a total volume of 10 μl for 3 h at 37°C according to the instructions of the manufacturer. Enzyme inactivation was carried out by incubation at 65°C for 20 min. The subsequent T-RFLP analysis was performed at the Genomics Technology Support Facility (http://genomics.msu.edu/ ; Michigan State University, East Lansing, Michigan). Briefly, the T-RFs were separated by capillary electrophoresis on an ABI Prism 3700 DNA analyzer. The DNA bands were automatically identified and sized using GeneScan software (Applied Biosystems, Foster City, CA) and comparison to internal lane standards. The relative abundances of individual T-RFs in a given pmoA PCR product were calculated based on the peak height of the individual T-RFs in relation to the total peak height of all T-RFs detected in the respective T-RFLP community fingerprint pattern. The peak heights were automatically quantified by the GeneScan software. To verify the assignments of T-RFs to our detected pmoA gene types, we also tested individual clones by T-RFLP analysis.

The T-RFLP results were highly reproducible. The coefficient of variation of the relative signal intensity of the T-RFs between different DNA isolations from the same soil sample ranged from 3 to 10.1%. The coefficient of variation of the relative signal intensity of the T-RFs between different PCRs from a single DNA sample ranged from 1 to 6.5%. Those variations are in the same range as those previously reported (28). The variations between different digests from the same PCR product and between different electrophoretic runs from the same digest were negligible. This is consistent with previous systematic evaluations of the T-RFLP method (38).

The relative abundance data were analyzed with a split-plot analysis of variance performed using the MIXED procedure in SAS (SAS Institute, Inc., Cary, NC). Means were estimated as least-square means, and the degrees of freedom were estimated using the Satterthwaite approximation. The data were arcsine-transformed before analysis.

The partial pmoA gene sequences determined in this study have been deposited in the EMBL, GenBank, and DDBJ nucleotide sequence databases under the accession numbers AY654669 through AY654732 .

Clone libraries were constructed using pmoA PCR products from three experimental plots: two with elevated CO₂, temperature, precipitation, and nitrogen (plots ID 5 and ID 60) and one with ambient levels of CO₂, temperature, precipitation, and nitrogen (plot ID 107). In total, 64 clones were analyzed (11 clones for ID 5, 21 clones for ID 60, and 32 clones for ID 107). Figure 1 shows the phylogenetic affiliation of all clone sequences analyzed in this study.

Five sequences formed a distinct clade (JR1) that was related to the “RA 14” clade, environmental sequence types that have been hypothesized to represent uncultured “high-affinity” methanotrophs capable of oxidizing methane at atmospheric concentrations (21, 25). The similarity in DNA sequence between JR1 and the “RA 14” clade was approximately 80%. Two sequences (“JR60-650-3” and “JR60-650-8”) grouped tightly with the “RA 14” clade.

Two further clades (JR2 with 22 sequences and JR3 with 14 sequences) were moderately related to each other but showed no close relationship to any cultivated methanotroph species. They grouped distantly to type I methanotrophs, with a DNA sequence similarity of approximately 72%. While representing distinct lineages, JR2 and JR3 branched together with the “WB5FH-A” clade, the second novel pmoA lineage suggested to represent atmospheric methane oxidizers (32). Although supported by a bootstrap value below 50%, the common branching point of JR2, JR3, and the “WB5FH-A” clade was supported by a tree calculated using the maximum-likelihood approach. In contrast, a neighbor-joining tree calculated from deduced amino acid sequences favored a common branching point of the “WB5FH-A” clade with type I methanotrophs; however, the bootstrap value was again below 50%. Phylogenetic analysis consistently suggests a common evolutionary origin for the sequence clusters JR2 and JR3 and the amoA gene of Nitrosococcus oceani (an ammonia-oxidizing bacterium that is capable of using methane as a carbon source and whose prevalence is thought to be restricted to aquatic systems).

Eleven pmoA sequence types (clade JR4) were closely related to Methylocystis parvus, a relatively well-characterized member of the type II methanotrophs. We also found six sequences (clade JR5) that grouped together with a novel pmoA lineage previously described as a diverged second pmoA gene copy present in various strains of type II methanotrophs (50). JR5 and the novel pmoA copy of one of the representative species (Methylocystis sp. strain SC2) had a DNA sequence similarity of 85%.

In summary, we identified pmoA gene types belonging to five different lineages within the phylogenetic radiation of the pmoA/amoA family. Three of these clades (clades JR1, JR2, and JR3) have DNA sequence similarities of 80% or less with previously described pmoA variants.

The ratio of nonsynonymous to synonymous nucleotide substitution rates (dN/dS) was determined for each novel clade. The overall dN/dS ratio (as calculated with the one-ratio model of codeml) was 0.11 for the type I side of the pmoA tree and 0.10 for the type II side (data not shown). The dN/dS ratios (as calculated with the freely varying model of codeml) along the branches leading to the three novel lineages (JR1, JR2, and JR3) were 0.10, 0.17, and 0.20, respectively (Fig. 2). The likelihood ratio test showed that these dN/dS ratios were not significantly different (P < 0.05) from the background dN/dS ratios in their respective sides of the tree (data not shown).

These dN/dS ratios could result from two possible scenarios. Clades JR1, JR2, and JR3 could have diverged recently, with insufficient time for their dN/dS ratios to reflect any change in functional state, or the three clades could have diverged earlier, with the low dN/dS ratios reflecting continued purifying selection. Using codeml, we estimated the relative divergence of the PmoA clades by estimating the likely numbers of synonymous changes along each of the three lineages, based on the assumption that synonymous changes are not acted upon by selection and accumulate steadily over time (55) (divergence times can also be roughly estimated by inspection of branch lengths in Fig. 1; however, these lengths reflect both synonymous and nonsynonymous changes). There were approximately 65 synonymous changes along the branch to JR1, 58 to JR2, 54 to JR3, and 109 on the branch leading to JR2 and 3. This is not substantially less than, for example, the 47 synonymous changes leading to Methylocapsa acidiphila, the 119 leading to the “WB5FH-A” clade, and the 55 leading to Type I methanotrophs (data not shown). Thus, the clades JR1, JR2, and JR3 do not appear to have diverged recently compared to other known pmoA clades, and so their estimated dN/dS ratios suggest that they are undergoing purifying selection, encoding functionally active proteins.

The pMMO and AMO genes are evolutionarily related (24), and at the amino acid level they share a number of highly conserved residues (25, 42, 52). Based on alignments of the predicted peptide sequences of the α subunits of 112 particulate methane monooxygenases (PmoAs) and 349 ammonia monooxygenases (AmoAs), Tukhvatullin et al. (52) identified residues common to both proteins. Ricke et al. (42) extended this analysis to include the second PmoA gene copy, PmoA2, present in many Type II methanotrophs (50). The inferred translation of the region amplified by the primers used in our study spans 16 of these highly conserved residues (Table 2). All members of novel clades JR2 and JR3 each had all 16 of these conserved residues. All members of JR1, JR4, and JR5 had 15 of the 16 conserved residues, with all but one member in each group also having the 16th residue. Among the residues common to both PmoA and AmoA, Tukhvatullin et al. (52) proposed a subset of seven that could potentially be the metal ligands of the active site. The translation of our amplified pmoA region spans three of these (residues E100, Y157, and H169), which are conserved in all Jasper Ridge sequences. A further set of four residues were identified as potential non-active-site metal ligands, which could additionally stabilize the peptide structure (52); our amplified region spans two of these (residues D182 and Y196), which are also conserved in all Jasper Ridge sequences.

In addition, Holmes et al. identified 21 residues that could distinguish PmoA from most AmoA sequences (25). Our pmoA amplicons spanned 16 of the putative PmoA/AmoA diagnostic residues. All of the clades we detected (with the exception of JR6) shared a high percentage of amino acid residues typical of PmoA (Table 2). JR4 had all 16 of the PmoA-specific residues, while JR5 and JR1 had 13 and 11, respectively, of the PmoA residues, and in all cases the mismatches were amino acids belonging to the same amino acid similarity group (22) as the conserved PmoA residue (Table 2). Both JR2 and JR3 shared 14 of the 16 PmoA residues. The two mismatches in JR2 and one in JR3 were in the same amino acid similarity groups as the conserved residues, while the other mismatch in JR3 was a perfect match in half of the sequences in this clade.

Figure 3 shows a representative community T-RFLP-profile and the assignment of the T-RFs to the sequence clusters detected in our study. All clones produced the T-RFs that were predicted based on the sequence information (data not shown). All pmoA clades determined by comparative sequence analysis could be consistently recovered by T-RFLP community analysis. JR2, JR3, and JR5 exhibited specific T-RFs (208 bp, 373 bp, and 349 bp, respectively), confirmed by in silico analysis of the publicly available pmoA gene sequences (combined with the sequences generated in this study). JR4 produced a T-RF of 245 bp as anticipated (this is the specific T-RF for the type II methanotrophs) (28). However, no specific T-RF could be generated for clade JR1 by use of MspI (i.e., the 80-bp T-RF generated by JR1 can also be produced by digestion of pmoA sequences from Methylococcus capsulatus and related species, as well as M. acidiphila). A T-RF of 34 bp was indicative for sequences belonging to the “RA 14” clade.

Although not confirmed by cloned sequences, our T-RFLP community profiles indicated the presence of various members of type I methanotrophs (e.g., T-RFs of 440 bp, 505 bp, and 511 bp, with the latter two representing undigested pmoA sequence types without the MspI recognition site) (28), although in low abundance (generally less than 4% of the total). We can think of at least two possible explanations for the absence of pmoA sequences related to type I methanotrophs in our clone libraries, namely: (i) low relative abundance of the type I methanotrophs combined with nonexhaustive clone sampling, and (ii) cloning biases against type I sequence types. Recently reported discrepancies between the community composition of pmoA clone libraries and pmoA-based T-RFLP analysis (28, 40) suggests that such biases can be present. Given this possibility, we did not attempt to determine the response of type I methanotrophs to simulated global change in our study.

We generated T-RFLP community profiles of methanotrophs from all replicate treatments of our multifactorial climate change experiment (8 replicates of 16 treatments, for a total of 128 soil samples). Simulated global change did not significantly alter the number of T-RFs present (the phylogenetic richness of the methanotroph community) or the magnitude of Shannon, Simpson, or Berger-Parker diversity indices (34) calculated from the T-RFLP data. However, the simulated global changes did alter community composition. The relative abundance of type II methanotrophs (clade JR4) significantly decreased under elevated precipitation (F_1,24 = 7.89; P = 0.0068) (Fig. 4) and elevated temperature (F_1,24 = 4.12; P = 0.0469) (Fig. 4). However, these effects were not additive; i.e., there was a significant antagonistic interaction between precipitation and temperature (F_1,24 = 8.31; P = 0.0055) (Fig. 4) such that the effect of both treatments together was less than that expected from their individual effects. In contrast, the relative abundance of the novel methanotroph clade JR2 responded to simulated global change very differently (Fig. 5). Elevated precipitation and temperature increased the relative abundance of this clade, and there was a significant antagonistic interaction between elevated precipitation and temperature (F_1,24 = 13.48; P = 0.0012) (Fig. 5) as well.

The aim of the present study was to explore the methanotrophic diversity of a Californian upland grassland and to assess whether a shift in the methanotrophic community structure in response to simulated global change was detectable. We assessed methanotroph diversity in this study using a cultivation-independent approach, with pmoA as a molecular marker. To date, most studies involving pmoA-based analysis of methanotrophic populations have used the primer system A189F-682R. These primers also amplify amoA, which encodes the homologous subunit of the ammonia monooxygenase in nitrifying bacteria. Reverse primers that discriminate against the amoA (e.g., mb661) and highly specific primers with intended target specificity for the “RA 14” clade (e.g., 650R) have been applied as alternative methods for studying methanotroph diversity. Bourne et al. (10) tested these three primer sets in various soils and found that one primer combination alone was not sufficient to explore methanotrophic diversity.

We tested five different primer combinations (three single-round PCR assays and two nested PCR assays) in order to determine their potential to detect a broad range of methanotrophs in our grassland soil. When primer set A189F-682R was used, clone libraries created from the single-round PCR amplicons showed a high representation of amoA inserts. When primer set A189F-mb661R or A189F-650R was used, the clone libraries contained a large number of nonspecific inserts. Nested PCR, however, using the A189F-682R primer set in the first round and either reverse primer mb661 or reverse primer 650R in the second round, generated consistently high yields of pmoA amplicons, even in some soils for which single-round PCRs produced little or no pmoA amplification. In fact, all analyzed clones derived from nested PCRs were “pmoA positive.” The reverse primers we used detected different components of the methanotroph community. Primer 650R detected the clades JR1 and sequences from the “RA 14” clade, clade JR3 could be detected only with reverse primer mb661, and clades JR2, JR4, and JR5 were detectable with both reverse primers. Therefore, we used both reverse primers together in a multiplex (i.e., in the same reaction), nested PCR approach for the T-RFLP community analysis. This enabled us to simultaneously recover a broad range of distinct pmoA clades in single electrophoretic profiles for each sample.

We discovered a remarkably high diversity of pmoA gene types in our study (Fig. 1), including those closely related to the pmoA of known members of the class Alphaproteobacteria as well as gene types distinct from known species forming hitherto undescribed pmoA lineages. Within type II methanotrophs, we found sequences closely related to M. parvus (clade JR4), as well as the recently characterized type II pmoA gene copy (50) of Methylocystis sp. (clade JR5). Interestingly, the relative abundance of the T-RFs was consistently higher for JR4 than for JR5 in our T-RFLP profiles (data not shown), which agrees with the findings of Tchawa Yimga et al. (50) that not all type II methanotrophs possess this additional gene copy. We also discovered the clade JR1, which forms a distinct subgroup of the “RA 14” clade, the clade that has been putatively identified as atmospheric methane consumers (21, 25). This finding considerably expands the known depth of the “RA 14” clade and demonstrates that methanotrophs possessing this gene type are not restricted to forest soils. We did not detect the other putative atmospheric methane consumers, the “WB5FH-A” clade (32), although we did discover two novel clades (JR2 and JR3) which are distantly related to the “WB5FH-A” clade.

There are several lines of evidence that suggest that the three novel pmoA clades we discovered (JR1, JR2, and JR3) encode functional monooxygenases, with a primary substrate of methane rather than ammonia. All three of the novel clades had dN/dS ratios well below 1 (Fig. 2), evidence for purifying selection (55, 56). If these genes were nonfunctional copies, a lack of selection would result in nonsynonymous changes occurring at the same rate as synonymous changes, pushing the overall dN/dS ratio towards 1; how closely it approached 1 would depend on the divergence time of these clades. The dN/dS of the branches leading to all three novel clades, however, are not statistically different from the “background” dN/dS in the rest of each respective half of the pmoA phylogeny. In addition, the high number of synonymous changes along the branches leading to the three novel clades suggests that they did not diverge recently (see Results above), and thus their low dN/dS ratios suggest that their encoded proteins are expressed and functional.

The conservation of functionally diagnostic amino acid residues provides further evidence for retained function in the novel clades and for their substrate specificity for methane rather than ammonia. The novel sequences contain a very high percentage of those amino acid residues conserved in both methane and ammonia monooxygenases (42, 52). These conserved residues include those proposed to bind metal ions within the active site and at secondary stabilization sites (42, 52), as well as a majority of the previously identified PmoA-specific residues (25). Among the mismatched residues, almost all are in the same amino acid similarity groups as the PmoA-specific residues. The novel Alphaproteobacteria clade JR1 has the lowest number of perfect matches to putatively PmoA-specific residues (11 of 16) (Table 2) and has several putatively AmoA-diagnostic residues. However, two of these “AmoA-like” residues are, in fact, shared by several other PmoA clades. Furthermore, JR1 robustly clusters in the Alphaproteobacteria, within which there are no known amoA-containing members. Thus, the total evidence suggests that JR1 likely binds methane rather than ammonia. The novel Gammaproteobacteria clades JR2 and JR3 did not contain any AmoA-diagnostic residues. However, this picture is complicated somewhat by the fact that the only known ammonia-oxidizing bacteria within the class Gammaproteobacteria, the N. oceani-like clade, also lack many of the AmoA-diagnostic residues, and they are the closest phylogenetic relatives of JR2 and JR3 (Fig. 1). However, based on protein and inferred-translation alignments, there appear to be six sites that distinguish the N. oceani-like AmoA from the Gammaproteobacteria PmoA (Table 2) and from the enzyme encoded by JR2 and JR3. At position 71, the N. oceani-like clade contains an AmoA-diagnostic residue present in no known PmoAs. This residue is not present in JR2 or JR3. At five other sites, the N. oceani-like clade contains conserved residues distinct from known PmoAs and AmoAs; two residues are at PmoA-/AmoA-diagnostic positions, and three others are at positions conserved in all other PmoAs and AmoAs examined (Table 2), strongly suggesting functional relevance. None of these residues is present in JR2 or JR3. Finally, hydrophobicity plots of the consensus protein sequence of JR2 and JR3 show four transmembrane domains at positions identical to those of the Gammaproteobacteria PmoA consensus; in contrast, the fourth domain of the consensus for N. oceani-like AmoA is shifted 12 residues towards the C terminus, exactly matching the position of the corresponding hydrophobic domain of AmoA found within the class Betaproteobacteria (data not shown). Together, these sequence analyses suggest strongly that JR2 and JR3 are more likely to preferentially bind methane than ammonia.

It has been suggested that feedback between methane flux and climate change may be due to changes in the structure of the methanotroph community (31, 51); however, it is unknown whether realistic global changes have the potential to alter the community structure of methanotrophs. We used T-RFLP analyses of pmoA to provide a molecular profile of the methanotroph community and to determine if shifts in community structure occurred in response to simulated global change. We observed shifts in the relative abundance of both type II methanotrophs and the novel methanotroph clade JR2.

Type II methanotrophs decreased in relative abundance in response to increased precipitation (under ambient temperature) (Fig. 4, compare the open bars). Previous studies have reported decreased methane oxidation rates under increased soil moisture (1, 11, 54), possibly due to limitations on the diffusive transport of methane through the soil gas phase when soil moisture is high (31, 46). It is reasonable that reduced oxidation rates could result in the reduced relative abundance that we observed here, although this was not directly tested in our study. We also observed a significant decrease in the relative abundance of type II methanotrophs in response to increased temperature (under ambient precipitation) (Fig. 4, compare the open and hatched bars on the left). Although the diffusion of methane can be altered by temperature (31), and rates of methane oxidation are known to vary with temperature, the effect we observed is unlikely to be caused by the direct effects of temperature on methane supply or oxidation. The change in soil temperature in our plots due to the temperature treatment is negligible. However, the temperature treatment in our experiment has been reported to significantly increase soil moisture at the time of year at which we sampled, due to effects on the plant community that alter water loss from plant transpiration in the spring (60). It is thus plausible that the decrease in relative abundance we observed with increased temperature is due ultimately to the same mechanism as the decrease we observed with increased precipitation: an increase in soil water content. Indeed, soil water content was significantly correlated with the relative abundance of type II methanotrophs (P = 0.0178), while other factors (ammonium, nitrate, plant biomass, net primary productivity) were not.

In addition, we observed a significant interaction between precipitation and temperature, such that the combined effect of increased precipitation and temperature on type II methanotrophs was less than that expected by their individual effects. It is unclear why this might be. It is not due to nonadditive effects of temperature and precipitation on soil water content; there was not a significant interaction between these two factors in regard to soil water content in our study (data not shown). One possible explanation is that the negative effects of soil moisture on methane diffusivity are ameliorated at higher water contents by an increase in the proportion of anoxic microsites in the soil, leading to a net increase in methanogenesis. This could result in the combined effects of temperature and precipitation being less than that expected by their individual effects, if when combined they raise the soil water content to a level where the proportion of anoxic microsites is increased. This hypothesis could be tested in future work by comparing the relative abundance of type II methanotrophs at our site across years that vary naturally in precipitation.

The relative abundance of the novel methanotroph clade JR2 also responded to elevated precipitation and temperature, although in a manner opposite of that of type II methanotrophs. The relative abundance of JR2 increased in response to elevated precipitation and temperature (Fig. 5), rather than decreased, as observed for type II methanotrophs. Why might JR2 have responded so differently from classical type II methanotrophs? If the methanotrophs in JR2 are atmospheric methane “specialists,” as suggested by their association (although distant) with the “WB5FH-A” clade, then they might be expected to out-compete other methanotrophs under low-methane conditions. Such conditions could be present under conditions of relatively high soil water content, such as those resulting from increased precipitation or temperature, which would reduce the diffusion of methane into the soil. We observed a significant antagonistic interaction between elevated precipitation and temperature, such that the combined effect on the relative abundance of JR2 was less than that expected from their individual effects. Although it is unclear why this might be, it is possible that it could be due to same mechanism suggested for type II methanotrophs: simultaneous increases in temperature and precipitation increase soil moisture such that methanogenesis increases, increasing the methane supply and reducing the competitive advantage of JR2. Again, this is a testable hypothesis.

Our observations of significant interactions among global change factors are consistent with previous studies of global change. For example, Shaw and colleagues observed that antagonistic interactions among global changes could alter plant biomass at our site (48). Furthermore, Horz et al. (26) observed that both the abundance and community structure of ammonia-oxidizing bacteria at our site were altered by antagonistic interactions among global change factors, including interactions between temperature and precipitation.

Our study expands our understanding of the diversity of naturally occurring pmoA gene types. The high number of novel pmoA clades we detected was possible because of our use of combinations of different PCR primers in a nested and multiplex manner. Since most other pmoA-based studies have relied on the use of only one primer set, it is plausible that the novel clades we observed are present in other environments as well but have been overlooked due to the primer set used.

Using this approach, we not only discovered novel pmoA clades, which evolutionary and sequence analyses suggest are functional, but we also observed that at least one such clade responded to simulated multifactorial global change in a very different manner than classic type II methanotrophs. To our knowledge, this is the first study that demonstrates that significant changes in the community structure of methanotrophs can occur in response to multifactorial global change. It is not yet known how widespread such responses are, how such responses may vary through time, or the relationship between such changes and ecosystem function. Nonetheless, our results demonstrate that methanotrophs can be altered by global changes and that multifactorial experimental approaches may be necessary to fully assess the complexity of these responses.