Hidden diversity and potential ecological function of phosphorus acquisition genes in widespread terrestrial bacteriophages

Abstract

Phosphorus (P) limitation of ecosystem processes is widespread in terrestrial habitats. While a few auxiliary metabolic genes (AMGs) in bacteriophages from aquatic habitats are reported to have the potential to enhance P-acquisition ability of their hosts, little is known about the diversity and potential ecological function of P-acquisition genes encoded by terrestrial bacteriophages. Here, we analyze 333 soil metagenomes from five terrestrial habitat types across China and identify 75 viral operational taxonomic units (vOTUs) that encode 105 P-acquisition AMGs. These AMGs span 17 distinct functional genes involved in four primary processes of microbial P-acquisition. Among them, over 60% (11/17) have not been reported previously. We experimentally verify in-vitro enzymatic activities of two pyrophosphatases and one alkaline phosphatase encoded by P-acquisition vOTUs. Thirty-six percent of the 75 P-acquisition vOTUs are detectable in a published global topsoil metagenome dataset. Further analyses reveal that, under certain circumstances, the identified P-acquisition AMGs have a greater influence on soil P availability and are more dominant in soil metatranscriptomes than their corresponding bacterial genes. Overall, our results reinforce the necessity of incorporating viral contributions into biogeochemical P cycling.

Fig. 1. Sampling sites and selected main characteristics of soil phages under investigation.

a Geographic illustration of the soil sampling sites of this study on the map of China. Detailed information of each sampling site is provided in Supplementary Data 1. The darker green color on the map is indicative of areas where samples were taken, and the lighter green color represents areas where no samples were collected. The map of China was obtained from http://geo.datav.aliyun.com/areas_v2/bound/100000_full.json and was visualized by the R package sf. b Sample-based accumulative curves of all viral operational taxonomic units (vOTUs) recovered from soil metagenomes of this study and those vOTUs encoding genes responsible for microbial P-acquisition (referred to as ‘P-acquisition vOTUs’, shown in the inset). The accumulative curves were generated by a custom R script (available on GitHub) utilizing the R packages foreach and doParallel. Data (n = 500) were presented as mean values ± standard deviations. Teal or light blue traces represent 500 iterations (sample order randomizations), and red or blue points are means. c, d The kingdom- and family-level classifications of the P-acquisition vOTUs according to geNomad (bar chart) and PhaGCN (pie chart), respectively. For better visualization, only the percentages of the three most dominant families to the total numbers of vOTUs that could be classified at the level of family are shown. Detailed information of the vOTU taxonomy is provided in Supplementary Data 8.

Fig. 2. P-acquisition auxiliary metabolism genes (AMGs) identified in this study.

a Numbers of the AMGs involved in inorganic P solubilization, organic P mineralization, P transportation, and P starvation response regulation, respectively. The numbers of each kind of AMGs recovered from individual habitat types are indicated by different colors. Detailed information on the AMGs is provided in Supplementary Data 7. b A graph illustrating the increase in gene kind after adding the AMGs identified in this study (orange) to those reported previously (gray; detailed information on them is provided in Supplementary Data 9). c Genome organization diagrams of four representative vOTUs identified in this study (i.e., vOTU54, vOTU22, vOTU47, and vOTU24) and three representative public isolated phage genomes (i.e., GCA_002606105.1, GCA_004016045.1, and GCA_002956935.1). Due to the large genome size of GCA_002956935.1, only 30 genes upstream and downstream of the pstBACS cluster are displayed. Predicted open reading frames are colored according to VIBRANT and DRAM-v annotation functions. VOG, virus orthologous groups. Detail information on these viral genomes is listed in Supplementary Data 5, 6, 11.

Fig. 3. Functional validation of three representative P-acquisition AMGs.

a, e, i show the computational protein models of two pyrophosphatases (i.e., PPa1 encoded by the AMG ppa1 and PPa54 encoded by the AMG ppa54) and one alkaline phosphatase (i.e., PhoD22 encoded by the AMG phoD22), respectively. Helices and sheets are colored in a rainbow scheme (from the N terminus in red to the C terminus in blue). Detailed information on individual computational protein models is provided in Supplementary Data 12. b, f, j show the activities of PPa1, PPa54, and PhoD22, respectively, with comparisons between each phage-encoded protein and its corresponsive controls. P1: the commercial pyrophosphatase (PPa, 1000 U mL⁻¹) of Escherichia coli was used as a positive control for both PPa1 and PPa54. N: the total protein from the recombinant E.coli cells transformed with an empty pET28a vector was used as a negative control for all three phage-encoded proteins. P2: the commercial recombinant E.coli alkaline phosphatase (1000 U mL⁻¹) was used as a positive control for PhoD22. The dots overlaying each bar represent the corresponding data points. c, g, k show the effects of pH on the activities of PPa1, PPa54, and PhoD22, respectively. d, h, l show the effects of temperature on the activities of PPa1, PPa54, and PhoD22, respectively. Data presented in (b–d, f–h, j–l) were mean values ± standard deviations from three independent experiments (i.e., n = 3). Relative activities of a given phage-encoded protein shown in individual panels (c, d, g, h, k, l) were calculated based on the highest activity of that protein reported within the corresponsive panel.

Fig. 4. The numbers and relative abundances of P-acquisition AMGs in individual sampling sites.

a The average numbers of all kinds of P-acquisition AMGs detected in the five habitat types. Horizontal lines represent the medians, while the boxes represent the interquartile ranges of the first and third quartiles. The vertical lines indicate the maximal and minimal values. The dots overlaying each bar represent the corresponding data points. Different letters on the top of the bars indicate significant differences between habitat types assessed with the two-sided Wilcoxon test, and the P value indicates the overall difference among all habitat types assessed with the Kruskal–Wallis test. b–f The relative abundances of the four categories of P-acquisition AMGs (illustrated with the bar charts, see scale values on the X-axis of the left-hand side of each panel) and the numbers of all kinds of P-acquisition AMGs (illustrated with white circles, see scale values on the X-axis of the right-hand side of each panel) detected in individual sampling sites. The four categories of P-acquisition AMGs (short for solubilization, mineralization, transporter, and regulator, respectively) are indicated by different colors. Sampling sites are first grouped as per their habitat types [Farmland (b), Forest (c), Grassland (d), Gobi desert (e), and Mine wasteland (f)], and then those within the same habitat type are arranged according to their latitudes (from south to north).

Fig. 5. Global distribution patterns of the P-acquisition vOTUs and AMGs identified in this study.

a Map showing the sampling sites of the global soil study (from a published global topsoil metagenome dataset) and the numbers of P-acquisition vOTUs of our study that were also detected in individual sampling sites of the global soil study. Circles represent the sampling sites and are colored based on habitat types. Circle sizes reflect the numbers of P-acquisition vOTUs detected in the corresponding sampling sites. Circles at the same coordinates are stacked according to their sizes, with the largest one at the bottom. *, the value in the bracket following a given habitat type represents the percentage of samples with P-acquisition vOTUs in that habitat type. The world map was generated by the function map_ data (“world”) in the R package ggplot2. b Histograms showing the total numbers of the P-acquisition vOTUs detected in individual habitat types of the global soil study. c Histograms showing the numbers of the global soil samples where the eight kinds of P-acquisition AMGs carried by the vOTUs shared by our study and the global soil study were detected. Habitat types are indicated by different colors. Detailed information is provided in Supplementary Data 15,16.

Fig. 6. Relative influences of selected factors on soil P availability evaluated by aggregated boosted tree (ABT) models.

a ABT results for the five habitat types as a whole. b–f ABT results for farmland, forest, grassland, Gobi desert, and mine wasteland, respectively. TP total phosphorous, MAT mean annual precipitation, Prokaryotic P-gene relative abundance of prokaryotic P-acquisition genes in the metagenome, Phage P-gene relative abundance of phage P-acquisition genes in the metagenome.

Fig. 7. P-acquisition vOTU-host linkages and gene transcription profiles of P-acquisition AMGs.

a Pairwise P-acquisition vOTU-host linkages identified in this study. vOTUs are shown on the right-hand side of the panel and colored according to the kinds of P-acquisition AMGs that they encoded, while the hosts are shown on the left-hand side of the panel and colored according to their taxonomy. The number of linkages between individual vOTUs and their hosts are proportional to the sizes of those lines linking them. Detailed information of the P-acquisition vOTU-host linkages is provided in Supplementary Data 17. b Histograms showing the numbers of kinds of P-acquisition AMGs transcribed in individual habitat types. The numbers of the transcribed genes belonging to the four categories of P-acquisition AMGs in each habitat type are also indicated by different colors. c The transcript ratios of phage:host gene pairs related to the four categories of P-acquisition AMGs in individual samples. Each circle in the panel represents a sample and is colored according to its habitat type. RPKM, reads per kilobase per million mapped reads. Detailed information on the public metatranscriptomes used to generate the results in (b) and (c) was provided in Supplementary Data 18, while the numbers of the samples where these P-acquisition AMGs were detected are shown in Supplementary Data 19.

Fig. 8. A schematic model of phage auxiliary metabolism in bacterial P- acquisition.

In low-P terrestrial ecosystems, P-acquisition phages can impact their hosts’ P-acquisition processes via four scenarios, which are indicated by red (Scenario 1), purple (Scenario 2), green (Scenario 3), and orange (Scenario 4) arrows, respectively. Positive and negative regulations of P-acquisition AMGs are indicated by solid and dashed arrows, respectively. Seven proteins (i.e., PPa, PhoD, PhoR, PhoU, PstA, PstB, and PstC) encoded by the P-acquisition AMGs detected in the investigated metatranscriptomes (listed in Supplementary Data 18) are shown; whilst two additional proteins (i.e., PhoB and PstS) are shown and marked by a question mark, as the transcripts of phage phoB and pstS were not detected in the metatranscriptomes. The PhoB with a capital letter ‘P’ indicates that it was phosphorylated. The substrates of PPa (i.e., pyrophosphate) and PhoD (i.e., phosphomonoester) are shown, respectively. For simplicity, only the outer and inner cell membranes for Gram-negative bacteria are shown. IM inner membrane, OM outer membrane, Pi phosphate.