Bacterial strains employed in this study are listed in Table 1 along with their characteristics relevant to methylotrophy. All have been previously isolated from Lake Washington sediment. The cultures were grown on solid media. The methanotrophs were cultivated on the nitrate minimal salts (NMS) medium in the atmosphere of methane (Dedysh and Dunfield, 2014; 25% methane, 75% air V/V). Methylotenera mobilis JLW8, Methylotenera sp. G11, and the Gram-positive methylotrophs were cultivated on the NMS medium supplemented with methylamine (0.2% W/V). The remaining methylotrophs were cultivated on the NMS medium supplemented with methanol (0.2% V/V). For the Rhodocyclaceae, vitamin B₁₂ was added at 0.5 μg L^-1 (Smalley et al., 2015). The non-methylotrophic heterotrophs were grown on diluted nutrient broth (BD Difco) medium (1/2 strength). Cultures were allowed to grow at room temperature (approximately 24°C) for 2–3 days, till they formed nice biofilms on the plates, after which the biomass was washed with the NMS medium, and optical density and/or cell counts were quantified, respectively, via spectrophotometry or flow cytometry (see Krause et al., 2017). The two exceptions were the Methylobacter strains, which were grown at 20°C, as they are psychrophilic (Kalyuzhnaya et al., 2015). The 50 strains were mixed at desired proportions (see below), and the master mixtures were incubated at 4°C for 48 h, to mimic starvation. Synthetic community 1 (SC1) master mixture was assembled by mixing the 50 strains in equal proportions, based on optical density (OD₆₀₀). The SC2 master mixture was assembled of the same strains, as part of a separate experiment, where they were mixed based on cell counts, as determined by flow cytometry. The main difference between SC1 and SC2 was in the proportion of the methanotrophs and the Bacillus strains, whose cells have larger size, thus they were less abundant in SC1 compared to SC2. Incubations were carried out at 18°C in 250 mL vials, with liquid-to-headspace ratio of 1–4, as previously described (Hernandez et al., 2015; Oshkin et al., 2015). Initial inoculum OD₆₀₀ was 0.1. Gas phase was refreshed daily according to the following scheme. (1) For the high methane high O₂ (HH) regimen, the headspace was flushed with air for 60 s (flow rate > 1 L/min), followed by removal of 50 mL of the headspace and replacement with 50 mL of methane; (2) for the high methane low O₂ (HL) regimen, the headspace was flushed with N₂ for 60 s (flow rate > 1.28 L/min), followed by the removal of 60 mL of headspace and replacement with 50 mL of methane and 10 ml of air; (3) for the low methane high O₂ (LH) regimen, the headspace was flushed with air for 60 s, followed by removal of 1 mL of headspace and replacement with 1 ml of methane; and (4) for the low methane low oxygen (LL) regimen, headspace was flushed with N₂ for 60 s, followed by removal of 2 mL of headspace and replacement with 1 mL of methane and 1 mL of air. When cultures reached OD₆₀₀ of approximately 0.5, they were transferred with 10-fold dilution into fresh medium. For the HH regimens, cells were transferred every 2–3 days, for the HL regimens, cells were transferred every 9–20 days, and for the LH and LL regimens, cells were transferred every 19–42 days. For the HH and HL regimens, cells from 11 time points were used for iTag profiling (Illumina sequencing of 16S rRNA gene fragment; Singer et al., 2016), for the LH and LL regimens, cells from five time points were used. For sample map and sampling schedule see Supplementary Table S1. Cells for DNA isolation (3 mL) were collected by centrifugation at 15,000 rpm for 5 min, and pellets were used immediately or stored at -20°C. Cells for RNA isolation (40 mL) were collected by centrifugation at 5,000 rpm for 30 min with the 1:10 volume of stop solution (5% buffer-saturated phenol in ethanol) added to cells prior to harvesting.

DNA was isolated using the FastDNA SPIN Kit for Soil (MP Biomedicals), according to the manufacturer’s protocol. Total RNA was isolated using the RNeasy minikit (Qiagen) as previously described (Chu et al., 2016). The purified RNA was tested for DNA contamination using 16S rRNA gene fragment PCR. Samples were stored at -80°C for subsequent analyses. iTag sequencing as well as transcript sequencing were performed at the DOE Joint Genome Institute (JGI) using standard JGI pipelines. Briefly, primers for hypervariable rRNA gene regions V4–V5 were used for iTag profiling (Parada et al., 2016). The primer sets contained Illumina linkers; a 12 bp barcode index; a pad region; a 0, 1, 2, or 3 base pair long spacer and the sequence-specific primer. Samples were pooled and run on two 2 by 300 bp Illumina MiSeq runs, essentially as described (Singer et al., 2016). RNAseq sequencing was carried out essentially as described (Beam et al., 2016).

iTag data have been archived with the IMG genome portal under JGI Project IDs 1105725 and 1105727. Transcript sequences have been archived with the IMG Genome Portal IMG genome IDs 3300009726, 3300009727, 3300009729, 3300009731, 3300009733, 3300009734, 3300009736, 3300009737, 3300009740, 3300009741, 3300009746, 3300009749, 3300009751, and 3300009752. Raw reads were aligned to the 50 respective reference genomes (see Table 1 for reference) using the Burrows–Wheeler alignment tool version 0.7.12-r1039 (Li and Durbin, 2009), using default parameters. Raw read counts are shown in Supplementary Table S2. Gene expression for organisms/genes recruiting large numbers of reads were normalized manually, by calculating numbers of reads per million of total reads mapping to a specific genome and normalizing per gene length (kb).

The UPARSE method was used for sequence processing and OTU clustering with USEARCH version 9.2.64 (Edgar, 2013). Clustering was performed at 97% sequence identity and chimeras were identified against the ChimeraSlayer reference database in the Broad Microbiome Utilities version r20110519 obtained from the UCHIME distribution (Edgar et al., 2011). Taxonomic assignments were made using the RDP Classifier from the Ribosomal Database Project version 2.12 (Wang et al., 2007). The samples were normalized so that the numbers of reads in each sample were equal.

Protein sequences were aligned using the CLUSTAL W algorithm, and phylogenetic trees were constructed using the maximum-likelihood method, as implemented in the MEGA7 software (Kumar et al., 2016). Statistical support was obtained from 1,000 bootstrap replicates (bootstrap values >50% are reported).

To address the questions posed, we first created SC1, by mixing 50 pure cultures of bacteria representing different functional guilds, all originating from Lake Washington sediment (Table 1). A single master mix was used to inoculate all vials. SC1 was incubated under three simulated environmental conditions: HH, HL, and LH (see section “Materials and Methods” and Supplementary Table S1 for details), in two replicates, nitrate being the nitrogen source in all cases. iTag sequence profiles were generated for 11 time points for each condition, reflecting community dynamics over, respectively, 27, 141, and 212 days (Supplementary Table S1). Under the HH condition, communities were mostly represented by two major OTUs, Methylophilaceae 1 (30–67% of total reads) and Methylomonas (12–62% of total reads), the proportion of Methylomonas increasing with time (Figure 1A). Methylosarcina OTU was detected in early samples (up to 23% of total reads), and the second Methylophilaceae OTU, Methylophilaceae 2 was detected across samples, its relative abundance declining with time (Figure 1A). An OTU representing Flavobacterium was detected in all samples, its relative abundance declining with time. To test for response to perturbation, a switch from HH to HL at day 28 was analyzed additionally. Community response included the switch between the Methylophilaceae OTUs, Methylophilaceae 2 gradually replacing Methylophilaceae 1, with overall increase in relative abundance of Methylophilaceae (Figure 1A).

The HL condition also selected for very simple community structure, dominated by Methylophilaceae 1 and Methylomonas, with minor presence of Methylosarcina, Methylophilaceae 2, and Flavobacterium (Figure 1A). An OTU representing Methyloversatilis was detected in some samples at over 2% abundance. Surprisingly, an OTU representing Acidovorax was detected in some samples at up to 5% abundance, even though Acidovorax was not intentionally added to the synthetic mix. Propagation of this organism must be a result of a minor contamination present in one of the frozen culture stocks.

Species composition was dramatically different under the LH condition, with Methylosarcina rapidly outcompeting all other methanotroph species, and with Methylophilaceae 2 outcompeting Methylophilaceae 1 (Figure 1A). An OTU representing Methylobacter was detected in SC1 samples only at very low abundance. Overall, the community dynamics and trajectories were very well replicated under each regimen (Figure 1A).

SC2 was built independently, of the same pure cultures revived from frozen stocks. With SC2, we expanded the range of simulated conditions to include nitrogen source as a variable, substituting ammonium for nitrate, and including an LL condition, a total of eight conditions. The HH and HL regimen communities were sampled at 11 time points, while the LH and LL regimen communities were sampled at five time points (Supplementary Table S1).

Community dynamics under the HH regimen were similar whether nitrate or ammonium were used as nitrogen sources (Figures 1B,C). In both cases, the Methylomonas types gradually outcompeted the Methylobacter types, and Methylophilaceae 2 gradually outcompeted Methylophilaceae 1 (Figures 1B,C). In contrast, community dynamics appeared to be different under the HL regimen, dependent on nitrogen source. While in both conditions, Methylomonas gradually outcompeted Methylobacter, in the nitrate-supplemented medium, Methylophilaceae were represented almost entirely by Methylophilaceae 1, and the sum of Methylococcaceae and Methylophilaceae made up to over 90 and up to 97% of the entire community (Figure 1B). However, in the ammonium-supplemented medium, Methylophilaceae were represented by Methylophilaceae 2, and the sum of Methylococcaceae and Methylophilaceae made up to no more than 80% of the community, other relatively abundant OTUs representing Flavobacterium, Methyloversatilis, and Acidovorax (Figure 1C).

Community dynamics under low methane (LH and LL conditions) were similar to each other, independently of O₂ concentration or nitrogen source. Methylomonas again outcompeted Methylobacter. However, in most cases, Methylophilaceae 1 prevailed over Methylophilaceae 2. An exception was one replicate under the LHNH₄ regimen, in which abundances of both Methylophilaceae OTUs increased with time. Under low methane, the sum of Methylococcaceae and Methylophilaceae made up approximately 80% of the total community, other abundant OTUs representing Flavobacterium, Methylosinus, and Acidovorax. Under the LHNO₃ regimen, we saw a significant population of Chryseobacterium (up to 13%), another organism that must have originated from a minor contamination of one of the frozen stocks. Notably, Methylosarcina was not detected in SC2 samples. Data for technical replicates, again, agreed very well (Figures 1B,C).

To summarize, while alternative major players within Methylococcaceae and Methylophilaceae were involved in community dynamics, especially dramatically contrasting the “high” versus “low” methane conditions, these two guilds were the dominant species in methane consumption, corroborating our previous results (Hernandez et al., 2015; Oshkin et al., 2015). As in our previous experiments, the Methylocystaceae species did not appear competitive under the “high” methane regimens. They were also gradually outcompeted by Methylococcaceae under “low” methane regimens.

Strain-resolved community compositions of SC2 microcosms were determined through transcript analysis of 14 samples (see Supplementary Table S1). These represented four time points for the two HH regimens (time series), and one time point for each of the remaining six regimens (Figure 2). Dynamics among the Methylococcaceae strains were dominated by Methylomonas strains under the HH regimens, with Methylomonas sp. LW13 gradually outcompeting Methylomonas sp. 11b and Methylomonas sp. MK1. However, the dynamics over time differed dependent on the source of nitrogen (Figure 2). Few reads were mapped to the Methylobacter genomes under these regimens, even though the Methylobacter OTU was relatively abundant at sampling points 1 and 3 in iTag analysis (Figures 1B,C). This suggests that Methylobacter species may have been transcriptionally inactive at the time of sampling. However, Methylobacter transcripts were detected under the HLNH₄ and the LH and the LL regimens, in proportions predicted via the iTag analysis. Methylosinus transcripts were relatively abundant under the LH regimens and also under the LLNO₃ regimen. Among the Methylophilaceae, Methylophilus species dominated over other Methylophilaceae under the HH regimens, Methylophilus methylotrophus Q8 gradually outcompeting M. methylotrophus 1. Methylophilaceae 7 was the dominant Methylophilaceae type under the HL and LH regimens when nitrate was supplied as the nitrogen source, and under both LL regimens. Methyloversatilis species outcompeted Methylophilaceae under the HLNH₄ regimen. Overall, abundance of the methanotroph transcripts as a fraction of total transcripts was higher in samples supplemented with ammonium (75–95%) compared to samples supplemented with nitrate (30–75%), likely due to transcription patterns responsive to the nitrogen source, as discussed below.

The most highly transcribed genes in the genomes of major species were determined through normalizing the number of transcript reads per gene length as a fraction of each species-specific transcriptome (Figure 3, Supplementary Figures S1–S4 and Table S2). Not surprisingly, genes encoding subunits of the particulate methane monooxygenase (pmoCAB) were most highly transcribed in the methane oxidizing species, both Methylococcaceae and Methylocystaceae, when nitrate served as the nitrogen source (Supplementary Table S2). Other methylotrophy functions were highly expressed in these species, as well as in the Methylophilaceae and the Rhodocyclaceae (Methyloversatilis). Remarkably, in all methylotrophs, some of the genes most responsive to specific experimental conditions were the genes encoding alternative MDH enzymes, MxaFI and XoxF. While the methylotroph genomes contain a single copy of each mxaF and mxaI (encoding the large and the small subunits of the Ca²⁺-dependent MDH; see Table 1 for genome reference) the number of xoxF genes (encoding the single subunit of the Ln³⁺-dependent MDH) is variable. The Methylococcaceae possess one (Kalyuzhnaya et al., 2015), the Methylocystaceae and Rhodocyclaceae possess two (Beck et al., 2015; Smalley et al., 2015), and the Methylophilaceae possess up to three (Beck et al., 2014). The phylogenetic relationships of MxaF and XoxF have been reviewed previously (Chistoserdova, 2011; Keltjens et al., 2014). The expression patterns were found very distinct dependent on whether nitrate or ammonium were supplied, or whether “low” or “high” concentrations of O₂ and methane were used (Figure 3 and Supplementary Figures S1–S4). Moreover, expression patterns differed among different species and among different MDH enzymes. While the Methylococcaceae, Methylocystaceae, Rhodocyclaceae, and the Methylophilus species expressed mxaFI genes at very high levels when “high” O₂, “high” methane, and nitrate were supplied, they expressed xoxF genes at higher levels when ammonium was supplied instead, or when methane was limited. Under the HLNO₃ regimen, genes for both enzymes were expressed at approximately equal levels. In contrast, Methylophilaceae other than Methylophilus (Methylophilaceae strains 7, 73s, 11) displayed higher expression of xoxF genes under all regimens. For example, at an early sampling point under the “high” methane, “high” O₂, nitrate, while Methylomonas, Methylobacter, Methylophilus, and Methyloversatilis were expressing genes for MxaFI at higher levels, non-Methylophilus Methylophilaceae species were expressing one or two genes encoding XoxF-type MDH enzymes at higher levels (Figure 4, Supplementary Figures S3, S4 and Table S2). Differential mxaF expression in response to the nitrogen source was confirmed on pure cultures of selected species, by real-time quantitative PCR (Supplementary Figure S5). It is noteworthy that the xoxF genes expressed by different community partners belong to two different phylogenetic clades (Chistoserdova, 2011), while Methylococcaceae, Methylocystaceae, and Rhodocyclaceae possess and express xoxF5, the Methylophilaceae possess and express xoxF4, and these genes are so far unique to this latter group (Figure 5).

In addition to the multiple genes for the MDH enzymes, genes for multiple cytochromes were highly and differentially expressed. Cytochrome MxaG, whose gene (mxaG) is cotranscribed with mxaFI genes, has been previously established as the electron acceptor from the MxaFI type MDH (Williams et al., 2006). While electron acceptors from XoxF enzymes have not been experimentally demonstrated, genes encoding cytochromes are colocalized on the chromosomes of alphaproteobacterial methylotrophs and Methylophilaceae (Beck et al., 2014, 2015). The arrangement is especially complex in Methylophilaceae, as they encode multiple cytochromes that belong to two phylogenetic lineages, XoxG1 and XoxG2, the former sharing homology with MxaG cytochromes, and the latter sharing homology with the cytochrome with a proposed function in methylamine oxidation (Kalyuzhnaya et al., 2008; Figure 6). In each case, genes for cytochromes representing the two phylogenetic groups were found co-expressed with specific xoxF genes, suggesting that they may be the natural electron acceptors from XoxF enzymes. Note that the coexpressing xoxF and xoxG genes did not always colocalize on the chromosomes (Figure 7).

While no cytochrome genes were present in the vicinity of the xoxF gene on the chromosomes of the Methylococcaceae, in each species we were able to identify a highly expressed gene, tentatively designated xoxG4, predicted to encode a novel cytochrome, with no sequence homology to either alphaproteobacterial xoxG3 genes or the two alternative Methylophilaceae genes, xoxG1 and xoxG2 (Figure 6). While XoxG4 is likely an electron acceptor from XoxF5 in the Methylococcaceae, xoxG4 transcription patterns differed from the ones of xoxF5 genes (Figure 3 and Supplementary Figure S1), suggesting that the role of this cytochrome is complex and may extend beyond methanol oxidation. Investigation of the properties of XoxG4 will be pursued in future studies.

Overall, for the Methylococcaceae, significant differences were noted between the transcriptomes of nitrate-supplemented versus ammonium-supplemented microcosms. In the latter, a series of hypothetical genes were expressed at very high levels, some even higher than the pMMO or the MDH genes. Interestingly, nitrogen fixation functions were expressed at high levels in ammonium-supplemented microcosms, compared to nitrate-supplemented microcosms (Figure 3 and Supplementary Figure S1). This type of regulation was somewhat unexpected. However, it may be potentially explained by faster oxygen consumption in the ammonium-supplemented microcosms, which are characterized by higher abundance of non-methanotroph species, and thus induction of nitrogen fixation genes may be in response to hypoxia. Future research is required to test this hypothesis.

Other highly expressed genes in the Methylococcaceae were the ones for pilin (pilA), proposed to be involved in extracellular electron transfer (Holmes et al., 2017; Lovley, 2017), and for type VI secretion functions, proposed to be involved in warfare among bacterial species (Alcoforado Diniz et al., 2015; Cianfanelli et al., 2016). In the Methylophilaceae, pilA genes were also highly expressed. Multiple porin genes and genes for other outer membrane proteins were also highly expressed (Figure 3 and Supplementary Figures S1–S4), and these may be involved in exchange of metabolites, signaling molecules, or other factors. The roles of these will also be addressed in future studies.