Abstract

   Non-methanotrophic bacteria such as methylotrophs often coexist with
   methane-oxidizing bacteria (methanotrophs) by cross-feeding on
   methane-derived carbon. Methanol has long been considered a major
   compound that mediates cross-feeding of methane-derived carbon. Despite
   the potential importance of cross-feeding in the global carbon cycle,
   only a few studies have actually explored metabolic responses of a
   bacteria when cross-feeding on a methanotroph. Recently, we isolated a
   novel facultative methylotroph, Methyloceanibacter caenitepidi Gela4,
   which grows syntrophically with the methanotroph, Methylocaldum marinum
   S8. To assess the potential metabolic pathways in M. caenitepidi Gela4
   co-cultured with M. marinum S8, we conducted genomic analyses of the
   two strains, as well as RNA-Seq and chemical analyses of M. caenitepidi
   Gela4, both in pure culture with methanol and in co-culture with
   methanotrophs. Genes involved in the serine pathway were downregulated
   in M. caenitepidi Gela4 under co-culture conditions, and methanol was
   below the detection limit (< 310 nM) in both pure culture of M. marinum
   S8 and co-culture. In contrast, genes involved in the tricarboxylic
   acid cycle, as well as acetyl-CoA synthetase, were upregulated in M.
   caenitepidi Gela4 under co-culture conditions. Notably, a pure culture
   of M. marinum S8 produced acetate (< 16 μM) during growth. These
   results suggested that an organic compound other than methanol,
   possibly acetate, might be the major carbon source for M. caenitepidi
   Gela4 cross-fed by M. marinum S8. Co-culture of M. caenitepidi Gela4
   and M. marinum S8 may represent a model system to further study
   methanol-independent cross-feeding from methanotrophs to
   non-methanotrophic bacteria.

Introduction

   Microbial methane oxidation plays an important role in the global
   methane cycle. Aerobic methane-oxidizing bacteria (methanotrophs) are
   key players in aerobic and micro-aerobic environments. The development
   of stable isotope probing has highlighted that methane-derived carbon
   is incorporated not only by methanotrophs but also by
   non-methanotrophic bacteria (methylotrophs or others) in diverse
   environments, and suggests the global importance of cross-feeding
   interactions in methane oxidation [[37]1, [38]2, [39]3, [40]4].
   Methanol is synthesized in the periplasm of methanotrophs, and
   therefore, may diffuse out of the cell [[41]5]. Hence, it is generally
   speculated that the co-existence of methanotrophs and
   non-methanotrophic bacteria mainly involves the cross-feeding of
   methanol excreted by methanotrophs to non-methanotrophic bacteria
   [[42]2]. Methanol-dependent cross-feeding between methanotrophs and
   non-methanotrophic bacteria has been studied since the 1970s [[43]6,
   [44]7]. Using isolates from lake sediments, Krause et al. [[45]8]
   demonstrated that methanol excreted by a methanotroph is indeed the
   sole carbon and energy source for the co-existing obligate methylotroph
   that can utilize methanol or methylamine exclusively.

   On the other hand, some methanotrophs are also known to excrete some
   organic acids [[46]9]. It is speculated that these organic acids are
   derived from dead cells. Recently, Kalyuzhnaya et al. [[47]10]
   demonstrated that Methylomicrobium alcaliphilum of Gammaproteobacteria
   release organic compounds, such as formate and acetate, via
   fermentation-based methanotrophy. This suggested that cross-feeding of
   methane-derived carbon from a methanotroph to non-methanotrophic
   bacteria may not be mediated solely by methanol and may involve more
   diverse processes. However, this has not been extensively explored yet.

   Recently, we obtained a methane-oxidizing enrichment culture from
   marine sediments of Kagoshima Bay, Japan and isolated a novel
   facultative methylotroph, Methyloceanibacter caenitepidi Gela4 [[48]11]
   and a novel gammaproteobacterial methanotroph, Methylocaldum marinum S8
   [[49]12]. M. caenitepidi Gela4 grows on a wide range of substrates and
   can grow syntrophically with M. marinum S8 using methane added as the
   sole carbon source. The co-culture of M. caenitepidi Gela4 and M.
   marinum S8 would be a suitable model system for addressing the nature
   of a facultative methylotroph cross-feeding on a methanotroph.

   In this study, we aimed to reveal the metabolic pathways employed by a
   marine facultative methylotroph cross-feeding on a methanotroph using
   genome analysis, RNA sequencing (RNA-Seq), and chemical analysis.

Materials and methods

Bacterial strains and culture conditions

   M. caenitepidi Gela4 was grown in NMS-SP medium [[50]12] with 1% (v/v)
   methanol added as a carbon source at 36°C. M. marinum S8 was grown in
   NMS-SP medium with a headspace of methane and air (20:80, v/v). For the
   co-culture of M. caenitepidi Gela4 and M. marinum S8, each culture was
   harvested by centrifugation (4,200 × g for 20 min), washed twice with
   NMS-SP medium and grown together in NMS-SP medium of 40 mL at 36°C in a
   200 mL glass bottle under methane/air (20:80, v/v). Culture growth was
   monitored by measuring the optical density at 660 nm. Cells were
   harvested by centrifugation at 4,200 × g for 20 min. For DNA analysis,
   the cells were then stored at −20°C. For RNA analysis, the supernatant
   was discarded after centrifugation, and 1.5 mL of RNAlater (Ambion
   Inc., Austin, TX, USA) was added. Re-suspended cells were stored in a
   refrigerator overnight. The cells were then harvested and stored at
   −80°C before RNA extraction, as described below.

Cell counts

   M. caenitepidi Gela4 and M. marinum S8 cells are morphologically
   distinct as visualized by a microscope [[51]11, [52]12]. Therefore,
   growth of M. caenitepidi Gela4 and M. marinum S8 in co-culture was
   monitored by cell count via microscopy. Cells were taken periodically
   and counted by a haemocytometer.

Genome sequencing, assembly, and annotation

   Chromosomal DNA was extracted using an ArchivePure DNA tissue kit (5
   Prime, Gaithersburg, MD, USA) according to the manufacturer’s
   instructions. The draft genome sequence of M. caenitepidi Gela4
   comprising of two scaffolds was generated previously [[53]11]. The
   remaining gap was closed by primer-walking, using primers complementary
   to the contig ends, at Takara Bio Inc. (Shiga, Japan). Whole-genome
   shotgun sequencing of M. marinum S8 was performed using PacBio RS II
   (Pacific Biosciences, Menlo Park, CA, USA) according to the
   manufacturer’s protocol. Genomes were assembled into a single contig
   using Sprai pipeline version 0.9.9.19
   ([54]http://zombie.cb.k.u-tokyo.ac.jp/sprai/) and manual curation. The
   contig was polished with Quiver [[55]13].

   The Microbial Genome Annotation Pipeline (MiGAP) was used to annotate
   sequences [[56]14]. Then, gene names were then assigned to genome
   annotations based on the best hits of blastp searches (version 2.2.30+;
   E-value<1e-5) [[57]15] against the SwissProt database (2017_2), with
   several genes annotated manually. The percentage amino acid identities
   between the predicted sequences determined in this study were
   calculated using the pairwise distances program in MEGA version 7
   [[58]16]. The Circos package [[59]17] was used to generate Circular
   plots. Genes were assigned to Kyoto Encyclopedia of Genes and Genomes
   (KEGG) pathway modules using metabolic and physiological potential
   evaluator (MAPLE) [[60]18]. Only modules with at least three assigned
   genes and qValue < 0.1 were retained, leading to 92 KEGG pathway
   modules for M. caenitepidi Gela4 ([61]S1 Table) and 100 modules for M.
   marinum S8 ([62]S2 Table).

   The in-house Julia scripts were used to define the duplicated genes in
   each species. The Bio.jl Julia package
   ([63]https://github.com/BioJulia/Bio.jl) was used to translate the
   annotated CDS in each genome into amino acid sequences. All-against-all
   blastp search was performed using the translated amino acid sequences
   and blast+ (version 2.2.30+) (E-value < 1e-5). A duplicated gene
   network, with nodes representing gene products and edges representing
   pairs of gene products, with E-value < 1e-5 and > 95% identity, was
   constructed. Genes with at least one edge in the duplicated gene
   network were defined as duplicated genes.

RNA extraction and sequencing

   RNA-Seq was used to compare the transcriptome profiles of M.
   caenitepidi Gela4 grown in pure culture on methanol with cells grown in
   co-culture with M. marinum S8 on methane-supplemented medium. Two
   replicate cultures at late exponential phase were prepared for each
   culture condition. The NucleoSpin RNA kit (Macherey-Nagel, Düren,
   Germany) was used to extract total RNA. Subsequent sample processing,
   including rRNA subtraction, library prep, and sequencing, was performed
   at Takara Bio Inc. Ribosomal RNA was removed from total RNA using the
   Ribo-Zero rRNA removal kit (Epicenter Biotechnologies, Madison, WI,
   USA). An Agilent 2100 Bioanalyzer was used to verify sample quality.
   RNA-Seq libraries were created using the Illumina TruSeq RNA sample
   prep kit version 2. All cDNA samples were individually barcoded and
   sequenced in the same lane on the HiSeq 2500 platform. The number of
   cDNA sequencing reads generated per sample varied between 22.0 and 29.9
   million.

Bioinformatics analysis of the transcriptome data

   For mapping of RNA-Seq reads, the genome sequences of the two species
   were combined into two chromosomes in a FASTA file. RNA-Seq data were
   mapped to the genome by STAR (version 2.5.2b) [[64]19] with the
   parameters “—outFilterMismatchNmax 3—outFilterMultimapNmax
   10—outSAMprimaryFlag AllBestScore—outSAMstrandField None—alignIntronMax
   3—alignIntronMin 4”. BAM files were processed using samtools 1.3
   [[65]20]. Specifically, uniquely mapped reads were selected with the
   parameter “-q 255” and primary alignment reads were selected with the
   parameter “-F 256”. The featureCounts option in the subread package
   (version 1.5.1) [[66]21] was used to count the number of reads mapped
   to each gene. For stain S8 data, the parameters “-M–fraction” were also
   set to calculate fractional counts for multi-mapping reads. Using
   simulated RNA-Seq data composed of the sequencing reads generated from
   the complete genome sequences of the two species, we confirmed that the
   aforementioned data processing strategy enabled unambiguous separation
   of the sequencing reads into two species ([67]S3 and [68]S4 Tables).
   Specifically, using custom Julia scripts, we generated simulated
   RNA-Seq data of each strain as FASTA files by extracting 100-bp length
   sequences of regions with at least 1-bp overlaps with genes. The
   sequences were prepared for each region from both forward and reverse
   strands. Reads were mapped to the combined genome of two strains and
   the primary alignment reads were selected as described above. Then,
   reads were counted for each CDS using the FeatureCounts with the
   parameters "-M—fraction". Differentially expressed genes were
   identified using edgeR (version 3.16.5) [[69]22] by quasi-likelihood
   test with genewise negative binomial generalized linear models and
   analysed by multivariate analysis. Only genes with ≥ 10 raw read counts
   in at least one sample were considered for the differentially expressed
   gene analysis. The reads per million mapped reads (RPM) value was
   calculated for each gene by dividing the read count by the total number
   of mapped reads. The reads per kb of transcript per million mapped
   reads (RPKM) value was then calculated for each gene, by dividing the
   RPM value by CDS length (unit = kb). The data for M. caenitepidi Gela4
   were based on unique mapping analysis, while the data for M. marinum S8
   were based on a primary alignment because of extensive gene duplication
   in the genome.

Chemical analysis

   Triplicate samples for chemical analysis were taken at the beginning of
   the incubation period and at early and late exponential phase. Samples
   of 2 mL were withdrawn using a syringe and filtered through a 0.22-μm
   pore size filter. Samples for methanol and acetaldehyde analysis were
   immediately sealed in a glass ampoule. Samples for organic acids and
   ammonium were stored in a plastic tube. All of the samples were stored
   at −20°C before analysis. Methanol and acetaldehyde concentrations were
   analysed by the Agilent 7697A Headspace Sampler coupled to the Agilent
   6890N Gas Chromatograph/MSD system with selected ion monitoring mode.
   One millilitre of the sample was added to a 5 mL headspace vial, heated
   at 80°C for 5 min, and 3 mL of headspace vapour was injected into a
   split injector of the gas chromatograph. The DB-WAX capillary column
   (60 m × 0.25 mm × 0.50 μm, Agilent J&W Scientific, CA, USA) was used as
   a GC column. Helium was the carrier gas at a flow rate of 1.0 mL/min.
   The detection limits for methanol and acetaldehyde were 310 and 570
   nmol/L, respectively. Concentrations of organic acids were analysed via
   HPLC (Shimadzu Corp., Japan) with a conductivity detector. Shim-pack
   SCR-102H (300 mm × 8 mm, Shimadzu model, Shimadzu Corp., Japan) was
   used as a column to separate organic acids eluted in 5 mmol/L of
   p-toluenesulfonic acid at 40°C. The detection limits were 1 μmol/L for
   formate and acetate and 10 μmol/L for propionate and butyrate.

Real-time quantitative reverse transcription PCR

   To support the RNA-Seq analysis results, expression levels of
   acetyl-CoA synthetase gene (acsA), hydroxypyruvate reductase gene
   (hpr), and aldehyde dehydrogenase gene (aldA) were quantified by
   real-time quantitative reverse transcription PCR assay. Co-cultures was
   incubated in triplicate and harvested at early and late exponential
   phases. M. caenitepidi Gela4 cells grown with methanol or sodium
   acetate (0.3%) were collected at late exponential phase and also
   analysed as reference. Total RNA was extracted as described above. For
   cDNA synthesis, 400–500 ng of RNA was reverse-transcribed using the
   PrimeScript RT reagent kit (Takara Bio, Japan) according to the
   manufacturer’s instruction. Primers specific for acsA, hpr and aldA
   genes of M. caenitepidi Gela4 were designed using Primer3 [[70]23]. The
   primers for acsA gene were GelaACS-F (5ʹ - GAACTCGACCTTTGCGTATCCG-3ʹ)
   and GelaACS-R (5ʹ - TGCAGTTCGCCGACACGTTAAG-3ʹ). The primers for hpr
   gene were GelaHPR-F (5ʹ - GATGCCGAACTTCATCGTCACG-3ʹ) and GelaHPR-R (5ʹ
   - GTTGCCCGCAACGAACGCATC-3ʹ). The primers for aldA gene were Gald-F (5ʹ
   - CGACACCGTTGCCTATCATTTC-3ʹ) and Gald-R (5ʹ - TCCAGCACGCCATCAAAATC-3ʹ).
   PCR amplification was not detected with these three primer sets when M.
   marinum S8 genomic DNA was used as a template. Real-time PCR was
   performed with TB Green Premix Ex Taq II (Takara Bio, Japan) using a
   Thermal Cycler Dice Real Time System II (Takara Bio, Japan). Each
   real-time PCR was carried out in a 15-μl total volume reactions that
   contained cDNA (1μl), Mastermix, forward and reverse primers (0.4 μM
   each) and nuclease-free sterile water. The thermal cycler program
   consisted of 40 cycles of denaturation (95°C for 5s) and annealing and
   extension (55°C for 30s) after an initial denaturation at 95°C for 30s.
   To construct calibration standard curves, dilution series of positive
   control DNA were used. For positive control DNA, gene fragments were
   amplified with the primers described above from M. caenitepidi Gela4
   genomic DNA. PCR products were then separated by gel electrophoresis
   and extracted with a NucleoSpin Gel and PCR Cleanup kit
   (Macherey-Nagel, Düren, Germany).

Results

Complete genome of M. caenitepidi Gela4 and possible metabolic pathways

   We generated the complete genome of M. caenitepidi Gela4 by closing the
   gaps from the previously published draft genome [[71]9]. The M.
   caenitepidi Gela4 genome was obtained at a sequencing coverage of 76×.
   It consisted of a single circular chromosome of 3,424,964 bp, with a
   G+C content of 62.8% ([72]Fig 1A). The genome was found to contain one
   ribosomal RNA operon (16S–23S–5S), 45 tRNA genes and 3,342 coding
   sequences (CDSs) ([73]S5 Table).

Fig 1.

   [74]Fig 1
   [75]Open in a new tab

   Circular diagrams of the (A) M. caenitepidi Gela4 and (B) M. marinum S8
   genomes. Circles (outside to inside) depict (1) contigs; (2) gene
   expression levels in pure culture, log10 (RPKM + 1) (only for M.
   caenitepidi Gela4); (3) gene expression levels in co-culture, log10
   (RPKM + 1); (4) coding genes on the leading strand (red); (5) coding
   genes on the lagging strand (blue); (6) 16S rRNA genes (green); (7)
   tRNA genes (blue); (8) transposase genes (grey); (9) duplicated genes
   (E-value<1e-5, % identity >95.0) (red); (10) the GC content (1-kb
   sliding window); and (11) the GC skew (1-kb sliding window).

   In M. caenitepidi Gela4, pathways for oxidization of methanol to
   formate were identified. The gene clusters responsible for methanol
   oxidation (mxaFJGIRSACKLD, GL4_0421–0431) and the synthesis of the
   cofactor pyrroloquinoline quinone (PQQ; pqqBCDE) were identified
   (GL4_3155, 2349–2351). In addition, four genes encoding a methanol
   dehydrogenase large subunit similar to xoxF were found (GL4_0437, 1130,
   1360, and 1361). The amino acid sequences of four XoxF proteins shared
   70%–94% identity with each other. All genes coding for enzymes involved
   in tetrahydromethanopterin (H[4]MPT)-mediated (fae, mtdB, mch and fhcD;
   GL4_1783 and 3256, 1788, 1786, 1794) and tetrahydrofolate
   (H[4]F)-mediated (mtdA, fch and ftfL; GL4_3376, 3375 and 3380)
   formaldehyde oxidation were identified. Genes for two formate
   dehydrogenases (FDH), NAD-dependent-FDH (fdhA1, B1, C1, D1, and E1;
   GL4_0089, 1367–1370), and tungsten-containing FDH (fdhA2 and B2;
   GL4_3220 and 3221) were also identified.

   Methylotrophs assimilate carbon from methanol by using either the
   ribulose monophosphate (RuMP) pathway or the serine cycle [[76]24]. In
   M. caenitepidi Gela4, all genes necessary for the serine pathway, i.e.
   serine hydroxymethyltransferase (glyA; GL4_1969), serine glyoxylate
   aminotransferase (sga; GL4_3378), hydroxypyruvate reductase (hpr;
   GL4_3377), d-glycerate 2-kinase (gckA; GL4_3381), enolase (eno;
   GL4_2130), phosphoenolpyruvate carboxylase (ppc; GL4_3372), malate
   dehydrogenase (mdh; GL4_0105), malate-CoA ligase (mtk; GL4_3373 and
   3374) and malyl-CoA lyase (mcl; GL4_3371) were identified. In contrast,
   genes for the RuMP pathway, i.e. 3-hexulose-6-P synthase (hps) and
   6-phospho-3-hexuloisomerase (phi), were not detected. The gene for
   isocitrate lyase was not found, but genes for the ethylmalonyl-CoA
   (EMC) pathway were identified as a glyoxylate regeneration pathway.
   These gene composition patterns indicated that M. caenitepidi Gela4
   assimilates carbon from methanol through the serine pathway. Moreover,
   genes for glycolysis and the tricarboxylic acid (TCA) cycle were also
   identified.

Complete genome of M. marinum S8 and potential metabolic pathways

   We generated the complete genome of M. marinum S8 obtained at 102 ×
   sequencing coverage using PacBio RS II. The genome was found to consist
   of a single circular chromosome of 6,088,916 bp, with a G + C content
   of 58.7% ([77]Fig 1B). The genome contains 48 tRNA genes and 5,611 CDSs
   ([78]Fig 1B and [79]S6 Table). Notably, it possessed a high number of
   transposase genes and duplicated genes ([80]Fig 1B); 9.5% (539) of
   genes were identified as duplicate when the threshold of amino acid
   sequence identity was set at >95% ([81]S1 Fig, [82]S7 Table). In
   addition, the genome contains two ribosomal operons (16S–23S–5S).

   Among the duplicated genes were key genes involved in methane
   oxidization. Two sets of particulate methane monooxygenase operons
   (pmoABC; sS8_1787–1789 and 3975–3977) were found; pmoA and pmoC genes
   were identical, whereas the predicted pmoB gene products showed 99.5%
   amino acid identity. In addition, the genome was found to contain three
   additional pmoC genes (sS8_2529, 3418 and 4876), whose predicted
   products shared a 60%–91% amino acid identity with the PmoC protein
   encoded by the pmo operon.

   Four genes predicted to encode the large subunit of methanol
   dehydrogenase were also found. One of them is located in a gene cluster
   encoding methanol dehydrogenase (MDH) and accessory proteins
   (mxaFJGIRSACKLD; sS8_1426–1436). The other three are similar to xoxF
   (sS8_0935, 1979 and 2523). Four genes for PQQ biosynthesis were found
   in the same cluster (pqqBCDE; sS8_2216–2219). All genes coding for
   enzymes involved in H[4]MPT-mediated and H[4]F-mediated formaldehyde
   oxidation were also identified. M. marinum S8 genome contains genes
   encoding the complete RuMP pathway, i.e. hps (sS8_1135 and 1148) and
   phi (sS8_1134) ([83]S6 and [84]S8 Tables). For the cleavage part of the
   RuMP pathway, genes for both the Entner-Doudoroff pathway (sS8_4013,
   4180, 4181, 4179, 1800) and Embden-Meyerhof-Parmas pathway (sS8_4382,
   1130) were identified. Among genes for the serine pathway, M. marinum
   S8 genome contains serine hydroxymethyltransferase (sS8_0137), serine
   glyoxylate aminotransferase (sS8_1474), hydroxypyruvate reductase
   (sS8_1475), D-glycerate-3-kinase (sS8_1476), enolase (sS8_1801), malate
   dehydrogenase (sS8_4317) and malyl-CoA lyase (sS8_2223) but lacks other
   genes. The M. marinum S8 genome contained all Calvin-Benson-Bassham
   cycle genes, such as ribulose-1,5- bisphosphate carboxylase (sS8_2432),
   phosphoglycerate kinase (sS8_4941), glyceraldehyde-3-phosphate
   dehydrogenase (sS8_1782), triosephosphate isomerase (sS8_0238),
   fructose-1,6-bisphosphate aldolase (sS8_1130), transketolase
   (sS8_1131), ribulose-phosphate 3-epimerase (sS8_ 2460) and
   phosphoribulokinase (sS8_4630). NAD-dependent FDH (fdhA1, B1, C1, D1
   and E1; sS8_0312–0316) and tungsten-containing FDH (fdhA2 and
   B2–sS8_2524 and 2525) were also found. The genome was also found to
   contain genes for the TCA cycle.

RNA-Seq of M. caenitepidi Gela4 in pure culture and co-culture with the
methanotroph

   To reveal the effects of syntrophic growth with a methanotroph on the
   gene expression profiles of M. caenitepidi Gela4, RNA-Seq was conducted
   on M. caenitepidi Gela4 cells grown on methanol in pure culture and to
   those grown in a co-culture with the methanotroph M. marinum S8 in a
   medium supplemented with methane ([85]S2 Fig). The specific growth
   rates of M. caenitepidi Gela4 in pure culture and the co-culture were
   0.03 and 0.02 (h^−1), respectively. [86]Fig 2 presents (A) growth of
   the co-culture and (B) species composition at three different growth
   phases. Percentage of M. marinum S8 in co-culture increased from 21% to
   77% during growth. The complete genome sequences of the two species
   enabled unambiguous separation of the sequenced reads between the two
   species in the co-culture samples ([87]S3 and [88]S4 Tables). When M.
   caenitepidi Gela4 was grown in co-culture with M. marinum S8, 24% of
   the uniquely mapped transcriptome data originated from M. caenitepidi
   Gela4 and 76% originated from M. marinum S8 on average. Pearson
   correlation coefficients of the expression levels between replicates
   were high in both pure culture (0.998) and co-culture (0.992) ([89]S3
   Fig), confirming high reproducibility of the RNA-Seq data. Full gene
   expression datasets for M. caenitepidi Gela4 and M. marinum S8 are
   presented in [90]S5 Table and [91]S6 Table, respectively.

Fig 2.

   [92]Fig 2
   [93]Open in a new tab

   (A) Growth of the co-culture and (B) species composition at early,
   middle and late exponential growth phase. The gray bar and white bars
   represent number of M. marinum S8 cells and M. caenitepidi Gela4 cells,
   respectively. Error bars represent standard deviation of triplicates.

   Expression of the central metabolism genes in M. caenitepidi Gela4 for
   both culture conditions, presented as RPKM are shown in [94]Table 1.
   The expression of central metabolism genes in M. marinum S8 in
   co-culture with M. caenitepidi Gela4, presented as RPKM, is shown in
   [95]S8 Table. [96]Fig 3 shows the central carbon metabolism of M.
   caenitepidi Gela4 and genes that were significantly up- or
   downregulated.

Table 1. Expression profiles of major genes of the central metabolism of M.
caenitepidi Gela4 cells grown on methanol in pure culture and in co-culture
with M. marinum S8; log2FC refers to log-fold changes in expression, and FDR
refers to false discovery rate.

     Predicted gene product Gene RPKM
   (pure culture) RPKM
   (co-culture) log2FC (co-culture/pure culture) FDR
   Methanol oxidation
   GL4_0421 Methanol dehydrogenase large subunit mxaF 42779 40524 -0.08
   0.308
   GL4_0422 Methanol dehydrogenase small subunit mxaJ 8017 8518 0.09 0.330
   GL4_0423 Cytochrome c-L mxaG 7201 6546 -0.14 0.074
   GL4_0424 Methanol dehydrogenase, small subunit mxaI 18826 21130 0.17
   0.041
   GL4_0437 Methanol dehydrogenase xoxF1 1117 1368 0.29 0.001
   GL4_1130 Methanol dehydrogenase xoxF2 80 98 0.29 0.001
   GL4_1360 Methanol dehydrogenase xoxF3 1123 547 -1.04 0.000
   GL4_1361 Methanol dehydrogenase xoxF4 1918 1115 -0.78 0.000
   Serine cycle
   GL4_1969 Serine hydroxymethyltransferase glyA 523 571 0.13 0.097
   GL4_3378 Serine-glyoxylate aminotransferase sga 320 162 -0.98 0.000
   GL4_2615 Serine-glyoxylate aminotransferase sga2 1023 467 -1.13 0.000
   GL4_3377 Hydroxypyruvate reductase hpr 568 235 -1.27 0.000
   GL4_3381 D-glycerate 2-kinase gckA 78 89 0.18 0.021
   GL4_2130 Enolase eno 176 201 0.20 0.021
   GL4_3372 Phosphoenolpyruvate carboxylase ppc 248 153 -0.69 0.000
   GL4_0105 Malate dehydrogenase mdh 272 417 0.61 0.000
   GL4_3373 Malate-CoA ligase alpha chain mtkB 705 405 -0.80 0.000
   GL4_3374 Malate-CoA ligase beta chain mtkA 684 435 -0.65 0.000
   GL4_3371 Malyl-CoA lyase mcl1 971 638 -0.61 0.000
   Acetyl-CoA synthetase and EMC pathway
   GL4_0367 Acetyl-coenzyme A synthetase acsA 98 1181 3.57 0.000
   GL4_0014 Acetyl-CoA acetyltransferase atoB 749 812 0.11 0.200
   GL4_0015 Acetoacetyl-CoA reductase phaB 331 418 0.34 0.001
   GL4_3328 3-hydroxybutyryl-CoA dehydrogenase hbdA 371 517 0.48 0.000
   GL4_3337 Crotonyl-CoA carboxylase/reductase ccr 869 1097 0.34 0.000
   GL4_3336 Methylsuccinyl-CoA dehydrogenase mcd 283 293 0.05 0.593
   GL4_2084 Methylmalonyl-CoA epimerase Mcee 173 165 -0.06 0.430
   GL4_3338 Ethylmalonyl-CoA mutase ecm 161 218 0.43 0.000
   GL4_0095 Mesaconyl-CoA hydratase mch 188 181 -0.05 0.367
   GL4_3371 malyl-CoA lyase mcl1 971 638 -0.61 0.000
   GL4_0468 L-malyl-CoA/beta-methylmalyl-CoA lyase mcl2 50 237 2.22 0.000
   TCA cycle
   GL4_2103 Citrate synthase citA 250 403 0.69 0.000
   GL4_0055 Aconitate hydratase acnA 294 153 -0.94 0.000
   GL4_1694 Isocitrate dehydrogenase icd 346 185 -0.90 0.000
   GL4_0108 2-oxoglutarate dehydrogenase E1 component sucA 386 652 0.75
   0.000
   GL4_0109 2-oxoglutarate dehydrogenase E2 component sucB 265 585 1.14
   0.000
   GL4_0107 Succinyl-CoA ligase alpha chain sucD 132 700 2.40 0.000
   GL4_0106 Succinyl-CoA ligase beta chain sucC 217 673 1.63 0.000
   GL4_0100 Succinate dehydrogenase flavoprotein subunit sdhA 252 217
   -0.22 0.004
   GL4_0101 Succinate dehydrogenase iron-sulfur protein sdhB 296 274 -0.11
   0.091
   GL4_0098 Succinate dehydrogenase cytochrome b-556 subunit sdhC 677 779
   0.20 0.009
   GL4_0099 Succinate dehydrogenase hydrophobic membrane anchor protein
   sdhD 140 136 -0.05 0.585
   GL4_1597 Fumarate hydratase fum 131 104 -0.33 0.000
   GL4_0105 Malate dehydrogenase mdh 272 417 0.61 0.000
   Ethanol oxidation
   GL4_0818 Alcohol dehydrogenase adhA 78 56 -0.49 0.000
   GL4_0813 Aldehyde dehydrogenase aldA 222 6354 4.83 0.000
   [97]Open in a new tab

Fig 3. Predicted central carbon metabolic pathway of M. caenitepidi Gela4.

   [98]Fig 3
   [99]Open in a new tab

   Red and blue arrows indicate genes that were upregulated (log2FC > 0.3,
   FDR < 0.05) and downregulated (log2FC < −0.3, FDR < 0.05) in
   co-culture.

   Differential gene expression analysis revealed that 689 genes of M.
   caenitepidi Gela4 were significantly upregulated [log2FC (fold-change)
   > 0.30, false discovery rate (FDR) < 0.05], and 764 genes were
   downregulated in co-culture when compared with pure culture (log2FC <
   −0.30, FDR < 0.05; [100]Fig 4). To identify pathway modules enriched
   for differentially expressed genes, genes in M. caenitepidi Gela4 were
   assigned to KEGG pathway modules using metabolic and physiological
   potential evaluator [[101]18], and gene set enrichment analysis (GSEA)
   [[102]25] was performed using the pathway module annotations. Despite
   the low number of genes assigned to KEGG pathway modules (656 of 3,342
   coding genes), 5 and 7 modules were significantly upregulated and
   downregulated, respectively ([103]Fig 5, FDR < 0.01). Among the
   upregulated pathway modules were ‘Citrate cycle, second carbon
   oxidation, 2-oxoglutarate = > oxaloacetate’, and ‘Citrate cycle (TCA
   cycle, Krebs cycle)’, whereas among the downregulated modules were
   ‘Formaldehyde assimilation, serine pathway’ and ‘Cytochrome c oxidase,
   cbb3-type’. We also performed the hypergeometric test for association
   of KEGG pathway modules to up- or downregulated genes, which yielded
   similar results ([104]S4 Fig, FDR < 0.05).

Fig 4. Differential expression of M. caenitepidi Gela4 genes in pure culture
and co-culture.

   [105]Fig 4
   [106]Open in a new tab

   Red and blue points represent genes whose expression was upregulated
   (log2FC > 0.30, FDR < 0.05) and downregulated (log2FC < −0.30, FDR <
   0.05), respectively.

Fig 5. The results of gene set enrichment analysis on KEGG modules.

   [107]Fig 5
   [108]Open in a new tab

   The modules with FDR < 0.01 are shown according to the enrichment
   scores.

Down-regulation of methanol assimilation pathways in M. caenitepidi Gela4 in
co-culture

   In M. marinum S8, genes involved in methane oxidization, including the
   methane monooxygenase (pmo operon) genes, methanol dehydrogenase
   (mxaFJ), hexulose-6-phosphate synthase, NAD-dependent formate
   dehydrogenase, and fae, showed high expression ([109]S8 Table). This
   suggested that methane oxidization occurs in M. marinum S8 as expected
   during the co-culture.

   In M. caenitepidi Gela4, genes encoding enzymes for the transformation
   of methanol to formaldehyde (mxa operon) and formaldehyde oxidation
   (fae) were among the most highly expressed in both pure culture and
   co-culture, and no significant changes in the expression levels of mxaF
   and mxaJ genes were observed (false discovery rate (FDR) > 0.05)
   ([110]Table 1). This implied that methanol is a primary carbon source
   of the methylotroph in co-culture with the methanotroph. However,
   surprisingly, methanol was below the detection limit in both pure
   culture of the methanotroph or co-culture ([111]Fig 6A). Considering
   that methanol dehydrogenase subunit 1 precursor is a major protein in
   both acetate-grown and methanol-grown Methylobacter extorquens AM1
   [[112]26], high expression of mxa operon genes and fae may not
   necessarily imply growth on methanol.

Fig 6.

   [113]Fig 6
   [114]Open in a new tab

   Concentrations of (A) methanol, (B) acetate, and (C) formate in culture
   supernatants taken at the beginning of the incubation, and early and
   exponential phase of M. caenitepidi Gela4, M. marinum S8, and
   co-culture. ND: not detected. Error bars represent standard deviation
   of triplicate experiments. Culture supernatants of M. caenitepidi Gela4
   at the beginning of the incubation and early exponential phase was not
   analysed due to high concentration (1% in the liquid phase at the
   beginning) for the analysis.

   In line with undetectable concentration of methanol, many genes in
   pathways involved in carbon assimilation through formaldehyde (the
   serine pathway (sga, hpr, ppc and mtk), and dissimilation (H[4]MPT
   pathway and H[4]F pathway) were significantly downregulated in M.
   caenitepidi Gela4 grown in co-culture compared to that in pure culture
   ([115]Fig 3). In addition, ‘Formaldehyde assimilation, serine pathway’
   was among the downregulated KEGG modules ([116]Fig 5). These results
   highlighted a downregulation of methanol utilization in M. caenitepidi
   Gela4 in co-culture compared to that in pure culture.

Up-regulation of TCA cycle in M. caenitepidi Gela4 in co-culture

   Next, we focussed on additional genes in central carbon metabolism
   pathways that were found to be upregulated in M. caenitepidi Gela4 in
   co-culture. The most highly upregulated gene was aldehyde dehydrogenase
   (aldA; GL4_0813; log2FC 4.83). The gene encoding the acetate oxidation
   enzyme, acetyl-CoA synthetase (acsA, GL4_0367), was also significantly
   upregulated (log2FC 3.57, [117]Table 1), which would lead to increased
   acetyl-CoA synthesis. The gene encoding the proton-translocating
   pyrophosphatase (hppa) (GL4_1976), which is involved in the cleavage of
   pyrophosphate, the by-product of acetyl-CoA synthetase, was also
   upregulated.

   Some genes in the EMC pathway, which converts acetyl-CoA to glyoxylate
   and succinyl-CoA, were upregulated in M. caenitepidi Gela4 in
   co-culture. Genes encoding acetoacetyl-CoA reductase (phaB; GL4_0015),
   3-hydroxybutyryl-CoA dehydrogenase (hbdA; GL4_3328), crotonyl-CoA
   carboxylase/reductase (ccr; GL4_3337), ethylmalonyl-CoA mutase (ecm;
   GL4_3338) and l-malyl-CoA/beta-methylmalyl-CoA lyase (mcl2; GL4_0468)
   [[118]27] were significantly upregulated ([119]Fig 3, [120]Table 1).

   In the TCA cycle, genes for citrate synthase (citA, GL4_2103), malate
   dehydrogenase (mdh, GL4_0105), succinyl-CoA synthethase (sucC and sucD;
   GL4_0106 and GL4_0107) and 2-oxoglutarate dehydrogenase (sucA and sucB;
   GL4_0108 and GL4_0109) were upregulated ([121]Fig 3, [122]Table 1). The
   TCA cycle was also found to be upregulated according to the pathway
   enrichment analysis ([123]Fig 5 and [124]S4 Fig).

   Importantly, many of the above genes that were differentially expressed
   in M. caenitepidi Gela4 were consistent with the expression profile of
   the facultative methylotroph Methylobacterium extorquens AM1 grown on
   acetate as compared with growth on methanol [[125]26]. Specifically,
   upregulation was observed in the genes involved in the TCA cycle.
   Furthermore, genes encoding the proton-translocating pyrophosphatase
   (hppa) (GL4_1976) and the NAD(P) transhydrogenase (GL4_3073–3075),
   which is used to balance redox equivalents in the cell, were
   upregulated in both M. extorquens AM1 grown on acetate and M.
   caenitepidi Gela4 in co-culture. PEP carboxykinase (GL4_0586), which is
   operating in an opposing fashion to the serine cycle PEP carboxylase,
   was also upregulated in both cases. In addition, downregulation of
   NAD-dependent FDH (GL4_0089, 1367 and 1368) and upregulation of
   tungsten-containing FDH (GL4_3220 and 3221) was observed in both cases.
   Collectively, these results suggested that carbon assimilation via the
   EMC pathway and energy generation via the TCA cycle also occur in M.
   caenitepidi Gela4 co-cultured with M. marinum S8, consistent with
   observations using ^13C steady-state metabolic flux analysis in M.
   extorquens AM1 grown on acetate [[126]26].

M. marinum S8 excretes acetate and formate

   To search for possible compounds cross-fed by the methanotroph to
   complement our RNA-Seq analyses, we measured the concentrations of
   organic acids in the culture supernatants. Acetate was produced at
   <15.6 μmol/L in the pure culture of M. marinum S8, while that in the
   co-culture and M. caenitepidi Gela4 was <12.4 and <1.9 μmol/L,
   respectively ([127]Fig 6B). Formate was detected at a concentration of
   <26.2, <8.4, and <12.4 μmol/L in culture supernatants of M. caenitepidi
   Gela4, M. marinum S8, and the co-culture ([128]Fig 6C). Acetaldehyde
   was below the detection limit in most of the culture supernatants, but
   a trace amount (590 nmol/L) was detected in one of the triplicate
   samples from the late exponential phase of M. marinum S8. Propionate
   was below the detection limit in all of the culture supernatants.

Quantitative analyses of acsA, hpr, and aldA genes

   RNA-Seq analysis suggested that M. caenitepidi Gela4 utilizes compounds
   other than methanol, possibly acetate, in co-culture. To further
   investigate and to monitor if the metabolic pathway changes with growth
   phase, real-time quantitative reverse transcription PCR was performed
   for acsA, hpr of the serine pathway and aldA genes of M. caenitepidi
   Gela4 in co-culture at early and late exponential growth phases. Pure
   cultures of M. caenitepidi Gela4 grown with methanol and acetate were
   also analysed as references. Student t-test revealed that the ratio of