Diversity and abundance of phosphonate biosynthetic genes in nature

Abstract

Phosphonates, molecules containing direct carbon-phosphorus bonds, compose a structurally diverse class of natural products with interesting and useful biological properties. Although their synthesis in protozoa was discovered more than 50 y ago, the extent and diversity of phosphonate production in nature remains poorly characterized. The rearrangement of phosphoenolpyruvate (PEP) to phosphonopyruvate, catalyzed by the enzyme PEP mutase (PepM), is shared by the vast majority of known phosphonate biosynthetic pathways. Thus, the pepM gene can be used as a molecular marker to examine the occurrence and abundance of phosphonate-producing organisms. Based on the presence of this gene, phosphonate biosynthesis is common in microbes, with ~5% of sequenced bacterial genomes and 7% of genome equivalents in metagenomic datasets carrying pepM homologs. Similarly, we detected the pepM gene in ~5% of random actinomycete isolates. The pepM-containing gene neighborhoods from 25 of these isolates were cloned, sequenced, and compared with those found in sequenced genomes. PEP mutase sequence conservation is strongly correlated with conservation of other nearby genes, suggesting that the diversity of phosphonate biosynthetic pathways can be predicted by examining PEP mutase diversity. We used this approach to estimate the range of phosphonate biosynthetic pathways in nature, revealing dozens of discrete groups in pepM amplicons from local soils, whereas hundreds were observed in metagenomic datasets. Collectively, our analyses show that phosphonate biosynthesis is both diverse and relatively common in nature, suggesting that the role of phosphonate molecules in the biosphere may be more important than is often recognized.

Fig. 1.

pepM gene abundance in GOS metagenomes and IMG/M microbiomes. (A) Boxplot of prokaryotic genome equivalents for pepM occurrence by habitat type in GOS. Single black lines represent the median value for environments that were sampled only once. No significant difference was found in the relative pepM abundance across GOS habitats (P = 0.9328, Kruskal Wallis test applied to habitats with more than one sampling site). (B) Boxplot of percentage of prokaryotic genome equivalents for pepM occurrence by ecosystem type in IMG/M microbiomes. Relative pepM abundance across various ecosystems differed significantly (P < 3.3 × 10⁻⁵, Kruskal Wallis test applied to categories with more than one sampling site). In A and B, the number of sampling sites for each type is shown in parentheses. (C) Distribution of predicted prokaryotic phyla for pepM homologs identified in GOS metagenomes. (D) Distribution of predicted prokaryotic phyla for pepM homologs identified in IMG/M microbiomes. Phyla that account for <1% are grouped in “Others.”

Fig. 2.

Analysis of PepM and phosphonate gene clusters from microbial genomes. (A) Maximum-likelihood tree of PEP mutase sequences from NCBI and 25 actinomycete strains from this study. The tree was calculated with the FastTree program (53) with default settings. The tree is rooted with 2-methylisocitrate lyase sequence (NP_286072). The branch is colored based on the source of pepM sequences: 25 actinomycete isolates from this study (red), NCBI archaeal genomes (purple), NCBI bacterial genomes (blue), and NCBI genomic fragments (green). Selected pepM groups are highlighted by color shading, with the number of sequences within a group shown in parentheses. Known phosphonate compounds are indicated and marked with an asterisk. The phosphonate biosynthetic loci for 25 actinomycete strains from this study are shown in SI Appendix, Fig. S2–S26. The phosphonate gene cluster for each organism is listed at http://file-server.igb.illinois.edu/˜xyu9/Dataset_S4._Phosphonate_gene_clusters.html. and is shown in the same order as in the tree. (B) Phosphonate gene cluster similarity as a function of PepM identity. Phosphonate gene cluster similarity was calculated using the fraction of homologous genes shared by two gene clusters. PepM identity was calculated using pairwise deletion of missing sites across the entire PepM alignment. Gene-cluster similarity measures were binned by PepM identity at intervals of 0.02 and plotted with a standard boxplot. The line shown is a linear regression over the PepM identity range 0.6–1.0, with equation ŷ = −0.86 + 1.6x; R-square = 0.74. The correlation has a P value of 2.2 × 10⁻¹⁶. The full dataset is plotted in SI Appendix, Fig. S28.

Fig. 3.

Rarefaction analysis of amino acid sequences of pepM identified from IMG microbial genomes, IMG/M microbiomes, GOS metagenomes, and soil pepM clone libraries. Rarefaction curves are shown for OTUs with differences not exceeding 16%.