A recently proposed re-organization of the taxonomy of the pink-pigmented facultative methylotrophs suggested renaming the extorquens clade with the genus name “Methylorubrum” (Green and Ardley, 2018). However, others have pointed out that the remaining Methylobacterium clades are not monophyletic, leading to the suggestion that the entire group be returned to the name Methylobacterium, until more data is gathered to support splitting the genus (Hördt et al., 2020; Alessa et al., 2021). Therefore, in this manuscript, we will refer to the formally described Methylorubrum species using that name, but to describe the larger monophyletic group that is the focus of this study we will use the term Methylobacterium.

The genetic potential of Methylobacterium species (Kelly et al., 2014) (Supplementary Table 1) for the degradation of aromatic compounds was assessed through analysis of genomes available on the IMG/M database (Chen et al., 2019). We surveyed all the genomes that were classified in the genera Methylobacterium or Methylorubrum as of April 2020. Of the 134 genomes returned, 3 were omitted from analysis: Methylobacterium sp. CCH7-A2 is classified by the Genome Taxonomy Database (Parks et al., 2018) as a Bosea sp. and by IMG as forming a clique with Porphyrobacter donghaensis; Methylobacterium sp. ZNC0032 is also annotated as a Bosea sp. and does not cluster with the Methylobacterium (included in Figure 1 and Supplementary Table 1 but not in further analyses); and Methylobacterium platani JCM 14648 has 2 genome entries in IMG, of which we used the earlier, more complete assembly. We should note that we found that a genome currently named Streptomyces purpurogeneiscleroticus NRRL B-2952 is identified as Methylobacterium by GTDB and in IMG forms a clique with several Methylobacterium species, including Methylobacterium fujisawaense, Magnaporthe oryzae, and Methylobacterium phyllosphaerae; however, it does not carry the vanA gene and was not included in this analysis. Finally, it has been noted that the published genome of Methylobacterium organophilum used in this study may contain contaminants and a resequencing places it in closer relation to the extorquens (B) clade than shown here (Alessa et al., 2021), but it is not one of the vanA-carrying organisms of interest. In total we analyzed 131 genomes, described in Supplementary Table 1.

Genes of interest were identified by Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology ID: vanA (vanillate monooxygenase alpha subunit, KO:K03862); pcaG (protocatechuate 3,4-dioxygenase alpha subunit, KO:K00448); pobA (p-hydroxybenzoate 3-monooxygenase, KO:K00481); and ech (trans-feruloyl-CoA hydratase/vanillin synthase, KO:K18383). To identify ech genes, we also searched for IMG Term Object ID 8865, vanillin synthase, and we used IMG’s Conserved Neighborhood search tool to identify an ech homolog in that was syntenic across several Methylobacterium genomes but only sometimes annotated as Vanillin Synthase. Locus tags of all genes are given in Supplementary Table 1. Vanillin synthase, vdh (KO:K21802), was not studied because we were not easily able to identify homologs in all organisms. Known lignin-degrading organisms were added to the phylogenies as reference organisms, and Prochlorococcus marinus NATL2A and archaeal species were used as outgroups (Supplementary Table 2).

For the phylogeny in Figure 1, which shows all putative aromatic-degrading Methylobacterium strains in the context of selected other Methylobacterium and known lignin-degrading organisms, four gene alignments were concatenated: 16S rRNA (the RNA component of the small ribosomal subunit), rpoB (DNA-directed RNA polymerase subunit beta, KO:K03043), atpD (ATP Synthase F1 complex beta subunit, KO:K02112), and recA (recombination protein RecA, KO:K03553), all of which have been shown previously to be informative for phylogenetic analysis in Methylobacterium and closely related organisms (Gaunt et al., 2001; Menna et al., 2009; Zhang et al., 2012). Locus tags are given in Supplementary Table 1. The clade structure generally agrees with previous Methylobacterium phylogenies based on 16S rRNA, rpoB, or core genes, though some analyses find clade A not to be monophyletic (Green and Ardley, 2018; Alessa et al., 2021; Leducq et al., 2022).

Two additional phylogenies were constructed as Supplementary Information: one containing all the Methylobacterium strains surveyed for this study, using the rpoB gene (Leducq et al., 2022) with B. diazoefficiens USDA 110 as an outgroup; and one showing the relationships of the Rieske [2Fe-2S] oxygenase homologs. The rpoB tree includes five of the Methylobacterium strains carrying aromatic catabolism genes (sp. MIMD6, sp. 174MFSha1-1, sp. AM11, sp. YR596, and sp. CCH7-A2) for which the full-length 16S rRNA gene sequence was not available and therefore are not included in Figure 1. For 3 of the Methylobacterium strains in our study (Methylobacterium sp. CG08_land_8_20_14_0_20_71_15, M. phyllosphaerae JCM 16408, and M. thiocyanatum JCM 10893), no rpoB sequences were available and aromatic catabolism genes were also not detected; these strains do not appear in any of the phylogenies but are listed in Supplementary Table 1. The homolog tree was constructed using all the genes annotated as COG 4638 in IMG in the surveyed Methylobacterium strains, as well as vanA sequences from other species for reference.

All gene sequences were aligned using MUSCLE v. 3.7 (Edgar, 2004) (RRID:SCR_011812) on the CIPRES Science Gateway (Miller et al., 2010), using the following parameters at all iterations: ClustalW as the sequence weighting scheme, UPGMB clustering, SUEFF = 0.1, pseudo rooting, and sum-of-pairs scoring for refinement. Phylogenetic trees were generated using RAxML v. 8.2.12 (Stamatakis, 2014) (RRID:SCR_006086) on CIPRES using the RAxML-HPC2 on XSEDE tool, with the following parameters: 25 rate categories, rapid bootstrapping with GTRCAT model for 100 iterations, random seed = 12345 for bootstrapping and parsimony, and outgroups as given in Supplementary Table 2 (outgroups are not shown in figures). All alignments are provided in Supplementary Files 1–10.

We used GeneHood v.0.15.0-0 (Ortega, 2020) to assess synteny in the gene neighborhoods surrounding the vanA and pcaH genes in Methylobacterium species. For each gene of interest, the database locus tag from the MIST 3.0 database (Gumerov et al., 2020) was provided and up to 15 genes upstream and downstream of the gene were plotted (in genomes that are not fully assembled, some scaffolds did not have 15 genes on either side of the target gene). Phylogenetic trees, gene annotations, and color-coding were added manually using Inkscape v. 0.91 (Inkscape Team, 2003) (RRID:SCR_014479).

GC content in the Methylobacterium sp. AMS5 genome was calculated in R v.4.0.2 (R Core Team, 2018) (RRID:SCR_001905) using the alphabetFrequency function from Biostrings v.2.46.0 (Pagès et al., 2019) (RRID:SCR_016949); for plotting, a sliding window of 5 kb was used for the full genome and 500 bp for the catabolic island region. Genome signature difference was calculated using a custom script in R, using the formula for the δ* difference (Karlin, 2001) with tetranucleotides rather than dinucleotides.

All organisms were grown on MPIPES medium (Delaney et al., 2013) with the addition of 1× Wolfe’s Vitamins (Atlas, 2010), 25 μM LaCl₃ [as has been found to facilitate growth in some Methylobacterium spp. (Skovran and Martinez-Gomez, 2015)], and carbon substrates as specified. All conditions contained 16 mM of carbon [16 mM methanol, 4 mM succinate, 2.67 mM glucose, 2.29 mM benzoate, 2.29 mM p-hydroxybenzoate (PHBA), 2.29 mM protocatechuic acid (PCA), 2 mM vanillate, 1.6 mM ferulate, or 0.4 g/L Kraft lignin]. While some soil microenvironments may experience higher substrate concentrations and some lower, 16 mM of carbon is in the range that has been commonly used for M. extorquens to ensure that final population size scales linearly with substrate concentration for both single- and multi-carbon substrates (Lee et al., 2009). Methylobacterium sp. 4–46 was also tested on methanol without LaCl₃ to assess the effect of lanthanum on methanol growth. For each assay, one colony was inoculated from a culture plate into 5 mL MPIPES with 3.5 mM disodium succinate and grown 24 h to obtain a stationary-phase culture, then diluted 1:64 (v/v) in the same medium and grown again until stationary. Using that inoculum, the growth experiment was conducted as follows: each replicate was diluted 1:64 into MPIPES and aliquoted into several wells of a Costar 48-well tissue culture-treated plate (product #3548, Corning Inc., Corning, NY, United States) for substrates other than methanol, or into Balch-type glass culture tubes for methanol, along with 16 mM carbon substrate. Total volume was 640 μL per well in culture plates or 5 mL per tube in Balch tubes. Each strain was tested in biological triplicate (3 separate colonies). Balch tubes were sealed with butyl rubber stoppers to prevent escape of methanol.

Plates were incubated at 30°C in a LPX44 Plate Hotel (LiCONiC, Woburn, MA, United States) with shaking at 250 RPM. Optical density was assessed using a Wallac 1420 Victor2 Microplate Reader (PerkinElmer, Waltham, MA, United States), reading OD₆₀₀ for 0.4 s at intervals of between 2 and 5 h. Culture tubes were incubated in a TC-7 tissue culture roller drum (New Brunswick Scientific, Edison, NJ, United States) at a speed setting of 7, and OD₆₀₀ was measured in the culture tubes using a Spectronic 200 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). Outputs from the Wallac 1420 software were collated using Curve Fitter (Delaney et al., 2013); data cleaning and plotting were then conducted in R.

Growth of Methylobacterium strains on vanillate was initiated with stationary-phase MPIPES-succinate cultures as described above; these were diluted 1:64 (v/v) into 5 mL MPIPES with 2 mM vanillate and grown 2 days until stationary phase to allow cultures to acclimate to the substrate. This inoculum was diluted 1:32 into culture flasks containing 20 mL of MPIPES with 2 mM vanillate and incubated, shaking, at 30°C. Cultures were sampled every 4 h: OD₆₀₀ was measured in a SmartSpec Plus Spectrophotometer (Bio-Rad Laboratories, Hercules, CA, United States), and the supernatant of 500 μL of centrifuged culture was used for the measurement of formaldehyde using the colorimetric assay described by Nash (1953). All strains were tested in biological triplicate. Although the active component in the Nash assay used has been found to react to several multi-carbon aldehydes, we have no reason to believe these compounds could be produced by vanillate metabolism, given what is known about the pathway; moreover, reactions with all compounds are slower than the reaction with formaldehyde (Compton and Purdy, 1980).

Formaldehyde production of non-methylotrophs during vanillate growth was assayed similarly. Stationary-phase vanillate-grown cultures of P. putida KT2440 and Rhodococcus jostii RHA1 [USDA Agricultural Research Service Culture Collection (NRRL)] were diluted 1:64 (v/v) into 5 mL MPIPES medium in culture tubes, with sampling every 2 h. Vanillate was provided at 4 mM; although this was greater than in the Methylobacterium experiment, we have found that vanillate concentration has little effect on the generation of formaldehyde in the medium (Lee et al., 2021).

To assess the abundance of Methylobacterium-encoded vanA genes in the environment, we searched all published metagenomes available through the IMG/M portal on December 20, 2017. The “Function Search” tool was used to identify all metagenomes containing genes annotated with KO:K03862 (vanA). To increase our chances of observing meaningful ecological patterns, we restricted our analysis to metagenomes containing >100 vanA genes. Genes on scaffolds with a phylogenetic lineage assignment of Methylobacterium were designated “Methylobacterium vanA.” and the Ecosystem Type field in the sample metadata was summarized into a “Sample Type” designation for analysis (e.g., Phyllosphere and Phylloplane were classified as “Leaf,” Rhizosphere and Rhizoplane as “Root,” etc.).

To provide a tree-based assessment of vanA diversity in addition to the USEARCH-based IMG/M phylogenetic lineage assignment (Huntemann et al., 2016), we downloaded the Methylobacterium vanA gene sequences and placed them in the reference vanA phylogeny (described above). pplacer v.1.1 (Matsen et al., 2010) (RRID:SCR_004737) was used with HMMer v.3.1b.2 (Eddy, 2007) (RRID:SCR_005305), RAxML v.8.2 (Stamatakis, 2014), and taxtastic v.0.8.8 for Python 3 (RRID:SCR_008394). The resulting tree was paired with metadata from sequence scaffolds in Phyloseq v.1.22.3 for R (McMurdie and Holmes, 2013) (RRID:SCR_013080) for the final phylogeny. Using pplacer, we found that only 182 out of 348 genes actually fell within the Methylobacterium clades, whereas 31 clustered more closely with Azospirillum halopraeferens, and the remaining 135 were distributed among more distantly related organisms (Supplementary Figure 1). This may be due in part to the overall complexity of the Rieske-type monooxygenase family (Özgen and Schmidt, 2019). We therefore took the conservative approach of focusing only on the genes identified by both IMG and pplacer as Methylobacterium for the remainder of our analyses.

A separate analysis was carried out to assess the relative abundances of all Methylobacterium reads within each metagenome in the JGI IMG database. The “Phylogenetic Distribution: Genomes v. Metagenomes” tool was used; this tool conducts a BLAST search to find reads in a set of metagenomes that have similarity to a query genome (JGI IMG, 2017). 18 query genomes were used: Methylobacterium radiotolerans JCM 2831, Methylobacterium pseudosasiacola BL36, M. extorquens PA1, Methylobacterium sp. AMS5, Methylobacterium sp. Leaf88, Methylobacterium gosipiicola Gh-105, Methylobacterium indicum SE3.6, M. indicum SE2.11, Methylobacterium aquaticum MA-22A, M. platani JCM 14648, Methylobacterium variabile 16961, Methylobacterium tarhaniae DSM 25844, M. aquaticum DSM 16371, Methylobacterium sp. WSM2598, Methylobacterium sp. 4–46, M. nodulans ORS 2060, Microvirga flocculans ATCC BAA-817, and Enterovirga rhinocerotis DSM 25903 (genome accessions in Supplementary Table 1). A total of 33,967 metagenomes were searched; this search set comprised all metagenomes published in IMG as of January 10, 2021 that were returned by the following 6 search terms:

Metagenomes with a published gene count of 0 were excluded from analysis, but otherwise no attempt was made to sample evenly across ecosystems, or to screen metagenomes by size, assembly status, etc., as the goal was to be as complete as possible in the search of available data. Due to computational constraints of IMG, metagenomes were searched in groups of 5,000 or fewer and 1 genome at a time, and results were combined. The output of each search consisted of a data table including identity (and metadata) of each metagenome sample where the genome was found, and the number of reads that hit to the genome at the 90% similarity level.

To investigate the possibility of aromatic catabolism among Methylobacterium strains, we first searched 131 published genomes for genes encoding four key enzymes known to be involved in the degradation pathway of methoxylated aromatic compounds: vanA (vanillate monooxygenase alpha subunit, KO:K03862); pcaG (protocatechuate 3,4-dioxygenase alpha subunit, KO:K00448); pobA (p-hydroxybenzoate 3-monooxygenase, KO:K00481); and ech (trans-feruloyl-CoA hydratase/vanillin synthase, KO:K18383; we searched also for IMG Term Object ID 8865, vanillin synthase, and used genome synteny; see Section “Material and Methods”). We found that all four genes are indeed present within the genus Methylobacterium. They are primarily limited to two closely related clades, the aquaticum clade (C1) and the nodulans clade (C2); furthermore, they are present in nearly all the genomes we surveyed in these clades (Figure 1, Supplementary Table 1, and Supplementary Figure 2). However, there are exceptions. Of the genomes searched, Methylobacterium sp. AMS5 in the extorquens clade (clade B) and M. pseudosasiacola in the radiotolerans clade (clade A) also carry some aromatic catabolism genes (Figure 1).

Two members of clade C2 (sp. 4–46 and sp. WSM2598) appear to lack vanA but carry ech (Figure 1). It is possible that the ech homologs annotated as trans-feruloyl-CoA hydratase/vanillin synthase in fact carry out a different reaction in these organisms (Lohans et al., 2017); however, they are located adjacent other genes annotated as vanillin dehydrogenase (see below), which catalyzes the conversion of vanillin to vanillate. This suggests that Methylobacterium sp. 4–46 and sp. WSM2598 may have previously had the capacity to utilize the product of these reactions but have since lost vanA. Among the closest sister taxa to Methylobacterium, E. rhinocerotis has none of the four genes of interest, and M. flocculans has only pobA and pcaH.

While the absence of a gene in a draft-status genome is not always an indication that the gene is not there, several of the genomes without aromatic catabolism genes, including M. sp. 4–46 and M. extorquens AM1 and PA1, are complete (Vuilleumier et al., 2009; Marx et al., 2012) (Supplementary Table 1). The substantial number of genomes we observed without any evidence of aromatic catabolism genes suggested a broad phylogenetic pattern: the capacity to metabolize several lignin-derived aromatic compounds was likely acquired by clade C early in its divergence within the Methylobacterium, but members of other clades have acquired the genes more recently by horizontal gene transfer. We sought further support for this evolutionary hypothesis by analyzing the phylogenies of the genes themselves, examining the genomic contexts of vanA and pcaG, and assaying the ability of several strains to grow on methoxylated aromatic compounds in the lab.

The phylogenies of pobA, pcaG, and vanA were largely congruent with the phylogeny of conserved housekeeping genes within the Methylobacterium but not between Methylobacterium and their closest relatives, supporting the hypothesis that the four genes are ancestral to clade C but were not necessarily acquired together (Figure 2). Phylogeny also supported our hypothesis of the function of the vanA gene in Methylobacterium species: although vanA homologs, members of the Rieske [2Fe-2S] oxygenases, are known to carry out a variety of reactions (Kweon et al., 2008; Özgen and Schmidt, 2019), and Methylobacterium species have several such genes, only the genes annotated as vanA cluster with genes that have been characterized as vanillate monooxygenases in known lignin-degrading bacteria (Supplementary Figure 3). The ech phylogeny was the most difficult to interpret, as the Methylobacterium separated into three clades (Figure 2). It is possible that not all gene homologs in the phylogeny code for enzymes specific to ferulate; further biochemical and genetic work will be necessary to be certain. Notably, Methylobacterium sp. AMS5 and M. pseudosasiacola were more closely related to each other than to any other species in both their pobA and pcaG sequences, despite their phylogenetic distance from each other at the genome level. And while the pobA sequences of those two strains fell within those of other Methylobacterium, their pcaG sequences clustered with the Sphingomonas, further support not only for the hypothesis that these two disparate species both acquired the beta-ketoadipate pathway via horizontal gene transfer, but that they may have acquired it from the same donor or that one served as the donor to the other.

We further examined the relationships of vanA and pcaG across species by comparing patterns of synteny within the genome neighborhood for each (Figure 3). In all strains with pcaG, the gene was located within a cluster of beta-ketoadipate pathway genes, and five of these genes were present in all strains in the same order. Strikingly, we found no genes annotated as pcaI and pcaJ (beta-ketoadipate:succinyl-CoA transferase) or pcaF (encoding beta-ketoadipyl-CoA thiolase) within these operons, which are necessary for the final steps of the beta-ketoadipate pathway (Harwood and Parales, 1996; Buchan et al., 2004), leaving it initially unclear whether and how these organisms might be able to metabolize vanillic acid and related compounds. However, almost all strains contained an additional gene between pcaG and pcaB, annotated as an uncharacterized conserved protein DUF849 (COG3246, beta-keto acid cleavage enzyme). In all Methylobacterium strains where vanA was found, it was accompanied by the other two genes required to confer the ability for the demethoxylation of vanillate: vanB (encoding the other subunit of vanillate monooxygenase) and vanR (an AraC family transcriptional regulator). Some strains also had genes encoding a Major Facilitator Superfamily transporter specific to protocatechuate (pcaK) or vanillate (vanK) (Figure 3).

The genomic contexts of pcaG and vanA (Figure 3), as well as pobA and ech (Supplementary Figure 4), exhibited patterns that agreed overall with the phylogenies described above. Closely related strains shared common genes in the regions near the genes of interest, with differences among the strains increasing with increasing distance from the genes of interest. Also consistent with phylogeny, Methylobacterium sp. AMS5 and M. pseudosasiacola shared no commonalities with members of clade C or with each other in the neighborhoods surrounding any of the genes of interest. In addition, these were the only two organisms in which pobA and pobR were located adjacent to the pca gene cluster; this contrasts with their arrangement in the other Methylobacterium and in most previously described organisms, where the pca genes are not co-located with any of the other genes of interest in this study (Supplementary Tables 1, 2). However, PCA, the reaction product of the enzyme encoded by pobA, is the substrate for the gene encoded by pcaGH; the proximity of the genes therefore suggests that if these two organisms gained the pca genes via horizontal transfer, the pob genes may have been transferred simultaneously as part of a catabolic island. We found further support for the catabolic island hypothesis in the proximity of vanA and ech to the pca genes in Methylobacterium sp. AMS5.

The co-localization of all four aromatic catabolism genes of interest in the genome of Methylobacterium sp. AMS5 prompted us to examine the region more closely. AMS5 was originally isolated in 2011 from the stem of a hypernodulating strain of soybean and, like other organisms in the group not formally described, is currently named as a Methylobacterium (Anda et al., 2011), but is very closely related to the model organism M. extorquens PA1 and therefore likely to belong to the Methylorubrum group (Figure 1). We aligned the two genomes, as well as that of another clade member, Methylorubrum zatmanii PSBB041, at the region where the aromatic catabolism genes are located (Figure 4). The comparison revealed that the region appears to be prone to the insertion and excision of mobile genetic elements.

At the locus where M. zatmanii carries a gene encoding a putative IS5-type transposase (possibly non-functional due to a frameshift with an internal stop codon), AMS5 contained an additional 22-kb region: a catabolic island containing genes that appear to confer the full ability to degrade ferulate, vanillate, vanillin, p-hydroxybenzoate, and protocatechuate, flanked by two different transposase genes from the MULE superfamily. On the 3′ side of this region, the two strains share 11 genes (14 kb), many related to nucleotide sugar metabolism, that are also not present in PA1 (Figure 4), or in any of the other Methylobacterium species surveyed in this study. The corresponding locus in the PA1 genome carries only a putative site-specific integrase/recombinase, and this gene is present only in members of the extorquens clade (clade B). The gene adjacent to the integrase encodes UDP-N-acetylglucosamine 2-epimerase (wecB); the sequence appears to be complete in M. zatmanii (1,137 bp) but truncated in PA1 (missing 912 bp from the 3′ end). In all three genomes, the region of variability lies immediately at the 3′ end of the genes for tRNA-Pro and tRNA-Arg. Due to the extensive genome rearrangement that is common within the Methylobacterium (Vuilleumier et al., 2009; Lee and Marx, 2012), we were unable to identify a syntenic region in species of clades A or C for comparison.

The presence of the truncated wecB gene in PA1 suggests that the region of sugar metabolism genes may originally have been present but were excised, whereas the arrangement of genes in the catabolic island suggest that the island was acquired into a M. zatmanii-like genomic background. The origins of the catabolic island itself are uncertain; as shown above (Figure 2), there appear to be diverse phylogenetic origins represented among the different genes within the pathway, and BLAST search for the entire region in the NCBI Nucleotide database found no other organisms outside of the Methylobacterium containing all genes in the same order. We could detect no genes in the region relating to replication or conjugation as might be expected in an integrative and conjugative element (ICE) (Wozniak and Waldor, 2010). Furthermore, no significant difference in tetranucleotide composition or in GC content was found between the inserted region (GC% = 71.4) and the full genome (GC% = 71.1) (Supplementary Figure 5), indicating that the catabolic island was likely either transferred from an unknown closely-related organism, or has been present in the genome long enough for amelioration (Lawrence and Ochman, 1997).

We next wanted to assess whether Methylobacterium strains carrying genes for the catabolism of methoxylated aromatic compounds could indeed use those compounds as growth substrates. We chose 8 species from the M. aquaticum/nodulans clades, shown in Figure 5 (and Supplementary Table 1). All of these strains carry all four genes of interest except Methylobacterium sp. 4–46, which lacks vanA and so was predicted not to grow on vanillate or ferulate (Figure 1). We conducted a series of growth experiments on defined mineral medium with substrates of each of the four gene products of interest, as well as several other substrates for context: methanol (a simple C₁ substrate), succinate (a known growth substrate for all the strains tested), glucose (used by only some Methylobacterium species), benzoic acid [another aromatic acid with a degradation pathway separate from the beta-ketoadipate pathway (Moreno and Rojo, 2008)], and Kraft lignin (a soluble form of polymeric lignin). For most strains, growth on PHBA, PCA, vanillate, and ferulate was as expected for all strains given their genome content: that is, organisms with pobA grew on PHBA, those with pcaG grew on PCA, those with vanA grew on vanillate, and those with ech grew on ferulate.

There were three exceptions: M. variabile did not grow on vanillate or ferulate, and M. nodulans did not grow on ferulate. Furthermore, M variabile growth on vanillate was only modest (Figure 5 and Supplementary Tables 3, 4). However, because a number of factors may lead to lack of growth, further tests in different culture conditions would be necessary to further explore these strains’ ability to grow on these substrates. For example, all strains grew on succinate, indicating that the cultures were viable and capable of growth on the medium provided; however, M. variabile and M. platani did not show substantial growth on methanol within 30 h, possibly indicating that conditions were not optimal for methylotrophic growth, or simply that a longer incubation time was needed. Growth on benzoate was observed for M. platani, M. tarhaniae, M. variabile, and M. nodulans but not for other species, and growth on Kraft lignin was negligible for all species in the 80 h tested.

We were especially interested to ask whether methylotrophs are different from non-methylotrophs in their ability to cope with the formaldehyde produced during vanillate metabolism. As mentioned above, it has been recognized by multiple researchers that the formaldehyde released by vanillate monooxygenase is a burden to organisms growing on vanillate, and we have previously observed formaldehyde accumulation during vanillate utilization by P. putida KT2440 (Lee et al., 2021). We therefore measured formaldehyde in the medium during vanillate metabolism by P. putida and Rhodococcus jostii RHA1, two well-studied lignin-degrading bacteria that, unlike Methylobacterium, do not utilize C₁ compounds for growth (Hibi et al., 2005; Ahmad et al., 2011; Sainsbury et al., 2013; Otani et al., 2014; Ravi et al., 2017). We assayed each strain on vanillate as well as on PCA as a control compound, as PCA is the direct product of vanillate demethylation and therefore involves the same metabolism except for the effect of formaldehyde. In both strains, when growth occurred on vanillate, formaldehyde accumulated in the medium concomitant with the increase in optical density, peaking at concentrations of 0.94 mM (Figure 6A and Supplementary Table 5). When cultures entered stationary phase (approximately 10 h for P. putida, 15 h for R. jostii), formaldehyde concentrations began to decrease, ultimately returning to below the detection limit, as might be expected from known formaldehyde detoxification pathways. No formaldehyde was detected at any time during growth on PCA. Moreover, growth on PCA was faster: stationary phase was reached at 7 h for P. putida and 10 h for R. jostii. These results suggest that while both P. putida and R. jostii can oxidize formaldehyde, degradation is slower than production.

We conducted a similar experiment on the five strains of Methylobacterium showing the greatest growth on vanillate in our laboratory conditions: M. aquaticum, M. nodulans, M. tarhaniae, Methylobacterium sp. AMS5, and M. aquaticum MA-22A. We incubated all cultures with vanillate for 33 h, long enough for AMS5 to show marked growth and all other species to consume all the substrate and reach stationary phase. No detectable formaldehyde was produced at any time (Figure 6B and Supplementary Table 5). These results suggest that Methylobacterium species are able to consume the formaldehyde generated from vanillate demethylation as rapidly as it is produced, likely via the same pathways by which the formaldehyde generated from methanol is used for energy generation and biosynthesis (Marx et al., 2005; Crowther et al., 2008).

Aside from the environments in which they were isolated, there exists scant information on the ecological niches of the aromatic-degrading Methylobacterium clades—for instance, whether they are more likely than the other Methylobacterium to inhabit ecosystems rich in lignin, such as soil or the rhizosphere. We therefore sought to learn more about the prevalence and abundance of vanA-carrying Methylobacterium species in the environment by querying the publicly available metagenome datasets on the JGI IMG/M database. The distinctiveness of the Methylobacterium genes within the phylogeny of known vanA sequences (Figure 2) makes it possible to deduce phylogeny from DNA sequence and vanA is present as a single copy in most genomes in which it is found.

Our study set comprised 1,651 metagenomes, with a total assembled gene count of 5.60 × 10⁹; we retrieved 317,816 scaffolds carrying vanA, of which 348 had Methylobacterium or Methylorubrum as their IMG Phylogenetic Lineage Assignment (a frequency of 0.11%) (Supplementary Tables 6–8). We then used pplacer (Matsen et al., 2010) (RRID:SCR_004737) for phylogenetic placement of those reads onto a reference tree of vanA genes from known Methylobacterium, and found that 182 of them actually clustered with the Methylobacterium; the others fell outside the Methylobacterium and were therefore omitted from further study (Supplementary Figure 1). Of those 182 Methylobacterium vanA, 114 (63%) clustered with the M. nodulans vanA gene; 31 genes (17%) with Methylobacterium sp. AMS5, and 37 genes (20%) fell into the M. aquaticum-like clade (Figure 7). Overall, they were most commonly found in samples associated with soil (43%) and roots (40%); only 10% were found in leaf samples.

To place these findings in context, we carried out a survey for Methylobacterium metagenome reads across all ecosystem types. We used the Phylogenetic Distribution of Genes in Metagenome tool (Markowitz et al., 2012; JGI IMG, 2017; Chen et al., 2019) to query all published metagenomes for reads with 90% BLAST similarity to each of 18 genomes from the Methylobacteriaceae. We observed a pattern in the samples that had BLAST hits to the query genomes suggesting ecological differences among the clades (Figure 8, Supplementary Figure 6, and Supplementary Tables 9–11). For example, the hits to query genomes from clades A and B were generally found with higher frequency in metagenomes from leaf samples than from soil (for Methylobacterium sp. Leaf88 and Methylobacterium gossipiicola, respectively, 4 and 10% of hits were from the soil, whereas 76 and 76% of were from leaf samples). In contrast, hits to query genomes from clade C2 were found primarily in soil and rarely in leaf samples (for M. nodulans, Methylobacterium sp. 4–46, and Methylobacterium sp. WSM2598, respectively, 74, 73, and 49% of hits were from the soil, whereas 2, 1, and 4% were from leaf samples). Clade C1 differed from the others in that there were fewer hits to these query genomes (between 8 and 71 hits per genome, in contrast to 96–487 hits for clades A–B and 104–523 hits for clade C2), and those hits more evenly distributed across all the sample types surveyed. These results point to ecological differences between the two clades of Methylobacterium capable of aromatic catabolism (clades C1 and C2). However, in general they agree with the survey of Methylobacterium vanA metagenome reads, suggesting that these two clades are more likely to be found in soil and root environments, and less likely to be found in the phyllosphere, than other clades of Methylobacterium.