Timing the evolution of phosphorus-cycling enzymes through geological time using phylogenomics

Domain: earth_science
Published: 2024-05-02
Editorial outcome: major_revision

Abstract
Phosphorus plays a crucial role in controlling biological productivity, but geological estimates of phosphate concentrations in the Precambrian ocean, during life’s origin and early evolution, vary over several orders of magnitude. While reduced phosphorus species may have served as alternative substrates to phosphate, their bioavailability on the early Earth remains unknown. Here, we reconstruct the phylogenomic record of life on Earth and find that phosphate transporting genes (pnas) evolved in the Paleoarchean (ca. 3.6-3.2 Ga) and are consistent with phosphate concentrations above modern levels (>3 µM). The first gene optimized for low phosphate levels (pstS; <1 µM) appeared around the same time or in the Mesoarchean depending on the reconstruction method. Most enzymatic pathways for metabolising reduced phosphorus emerged and expanded across the tree of life later. This includes phosphonate-catabolising CP-lyases, phosphite-oxidising pathways and hypophosphite-oxidising pathways. CP-lyases are particularly abundant in dissolved phosphate concentrations below 0.1 µM. Our results thus indicate at least local regions of declining phosphate levels through the Archean, possibly linked to phosphate-scavenging Fe(III), which may have limited productivity. However, reduced phosphorus species did not become widely used until after the Paleoproterozoic Great Oxidation Event (2.3 Ga), possibly linked to expansion of the biosphere at that time.

Introduction
Phosphorus is one of the fundamental building blocks of life. It is an essential constituent of DNA, cell membranes, and ATP, but is quickly depleted by plants and microbes in marine and terrestrial settings worldwide. As a result, restricted phosphorus availability limits primary productivity in large areas of the continents and oceans1,2.

It is conceivable that Precambrian microbial communities were similarly limited by phosphorus, and if so, that may have constrained the growth and evolution of the biosphere3–6. However, geochemical estimates of the Precambrian marine phosphorus reservoir vary by several orders of magnitude, from less than 0.12 μM to 4,100 μM in the Archean and between 0.03 and 2.5 μM in the Proterozoic oceans7–16. It has also been proposed that local microscale environments could have achieved up to 100 mM orthophosphate to drive prebiotic chemical pathways17. These disparities arise from differences in sampling strategies and assumptions in computational models. Therefore, the extent to which phosphorus was available to microbes in the Archean and Proterozoic remains a topic of debate, especially regarding surface waters with fewer direct constraints.

Modern microbial communities cope with limited phosphorus availability using a rich variety of enzymes to source phosphorus from the environment, including a range of dissolved organic phosphorus (DOP) compounds and inorganic phosphorus species, collectively referred to as total dissolved phosphorus (TDP) (Fig. 1). These include orthophosphate and phosphate esters (where P is in a fully oxidised state; P(V)) as well as more reduced compounds in the form of phosphite (P(III)), phosphonates (likely P(III)) and hypophosphite (P(I))18–26. Phosphite and hypophosphite constitute up to 26% of DOP in some habitats24 and geochemical estimates suggest that phosphite in particular was prominent in the Paleoarchean ocean due to elevated input from meteorites and Earth’s more reducing redox state27,28. Furthermore, phosphite is a source of phosphorus for multiple microbial species22, and in some cases, it can serve as an electron donor for chemolithotrophs25,29,30. In contrast, phosphonates make up approximately a quarter of all high-molecular weight DOP in the global oceans31. They contain carbon chains and characteristically strong carbon-phosphorus (C-P) bonds which are not present in other reduced phosphorus molecules. They can be assimilated as a nutrient resource by a diverse range of bacterial and archaeal taxa across the Tree of Life and in several different habitats18,19,21–23. Importantly, genomic studies have found that phosphonate-catabolising pathways catalysed by CP-lyases, phosphonatases, and some dioxygenases are more common in phosphate-depleted waters20,32 and soils33, suggesting that they may represent a mechanism for coping with phosphate scarcity. This observation implies that the phylogenetic history of these enzymes can be used as a tool for reconstructing phosphate availability in the past.

Here, we perform phylogenetic analyses to investigate which phosphorus compounds were available for microbial communities over geologic time. More specifically, we reconstruct the evolutionary history of microbial genes involved in utilising phosphate, phosphonates, phosphites, and hypophosphites. To do this, we reconcile a time-calibrated tree of life with phylogenetic trees of phosphorus-utilising genes to estimate when these genes were gained, lost, duplicated, and vertically inherited (see “Methods”). We then infer past phosphate concentrations by comparing presence/absence and abundance patterns for each gene with phosphate concentrations in the modern ocean. Our results allow us to shed new light on the evolution of the biogeochemical phosphorus cycle and provide the first constraints on the availability of reduced phosphorus substrates.

Results
Tree of Life
To determine how ancestral microbes used the various phosphorus species considered in this study, we began by reconstructing a tree of life from 865 genomes, chosen to represent the full diversity of life and genomes harbouring ptxD. Our results (Fig. 2) are broadly consistent with previous studies. For example, Actinobacteria are sister to Firmicutes (ultrafast bootstrap support (UBS)=95) and Proteobacteria probably share a more recent common ancestor with Bacteroidetes (UBS=66) than with Elusimicrobiota (UBS=37, Fig. 2) although support values for the latter are low indicating a high degree of uncertainty. A similar relationship has been predicted by Martinez-Gutierrez and Aylward34 and Coleman, et al.35, but not Moody, et al.36.

To estimate when key groups of organisms evolved, we implemented one uncorrelated (UGAM) and two autocorrelated (LN and CIR) Bayesian molecular clock models on alignments of 16 ribosomal proteins that have been used previously to reconstruct the tree of life30. The resulting time-calibrated phylogenies made with CIR, LN and UGAM clock models suggest that the last universal common ancestor (LUCA) radiated into Bacteria and Archaea at ca. 4.35 Ga (confidence intervals range from 4.40–4.22 Ga), 4.38 Ga (confidence intervals range from 4.40–4.34 Ga) or 4.37 Ga (confidence intervals range from 4.40–4.29 Ga), respectively. These dates fall shortly after the first evidence for liquid water37,38, which was applied as an upper constraint for the origin of life.

All three molecular clock models also predict that crown Bacteria emerged at least 200 million years before crown Archaea (Table 1), consistent with previous studies36. However, if methanogenesis evolved earlier ≥3.46 Ga instead of ≥2.7 Ga, in line with new geochemical evidence39–41, crown Bacteria may not have pre-dated Archaea because the estimated origin of crown Archaea is pushed back ~0.5 billion years to 3.74 Ga (confidence intervals span 3.83–3.66 Ga) whilst the estimated origin of crown Bacteria remains unaltered at 3.80 Ga with LN clock model (confidence intervals span 3.92–3.68 Ga).

Distribution of phosphorus-cycling genes in the Tree of Life
We surveyed 865 genomes for 13 genes involved in the utilisation and production of reduced phosphorus molecules and phosphate. The results reveal a broad distribution of genes encoding phosphate sodium symporters (namely pnas) and active phosphate transport mechanisms (indicated by pstS) throughout 257 and 174 strains of our dataset respectively (Fig. 2). These encompass all major bacterial and archaeal phyla except DPANN and Patescibacteria (formerly CPR, Fig. 2). In contrast, genes for low-affinity phosphate import of orthophosphate ions bound to divalent metal ions, such as Mg2+, Ca2+, Co2+, Mn2+ and Zn2+ via pitA and pitH are absent or rare respectively. Only one of the 865 genomes harbours a pitH homolog (Fig. 2).

Reduced phosphorous metabolisms are present in a variety of bacteria, but most archaea cannot metabolise phosphonates, phosphite or hypophosphite (Fig. 2). Of these, bacteria from the Proteobacterial superphylum host a particularly broad variety of phosphonate-cycling genes, as has been documented in previous studies20,21. We also find that phosphonate production is rarer than its catabolism, (consistent with previous work18), because only 30 genomes harbour pepM and ppd to produce phosphonates from phosphate esters, whereas 94 genomes have at least one of phnX, phnA, phnZ, phnJ or phnM to catabolise phosphonates. Of all the phosphonate catabolism genes, phnJ and phnM encoding the broad-specificity CP-lyases are most common because they co-occur in 44 genomes compared to the substrate-specific phosphonate catabolism genes, phnX, phnW, phnA and phnZ which are each present in just 28, 25, 12 and 40 genomes respectively (Supplementary Table 1).

The distribution of these genomes in the Tree of Life (Fig. 2) reveals a broader taxonomic distribution of phosphonate-producing compared to phosphonate-catabolising organisms, which is consistent with previous work20.

Hypophosphite is more reduced than phosphonates and phosphite. As a result, only 4 genomes have htxB encoding the binding protein required for active hypophosphite import from the environment. Two of these genomes also have htxA for oxidising hypophosphite to phosphite, but a further 15 genomes have htxA without the binding protein for hypophosphite import (Fig. 2). these htxA genes may encode enzymes that oxidise formate instead, as biochemical analyses have found comparable rates of activity on these small organic molecules42.

Timing of gene events
To estimate when these phosphorus-cycling genes emerged and radiated across the tree of life, our molecular clocks were reconciled with Bayesian evolutionary trees of each phosphorus-cycling gene to ascertain the timing of speciation, duplication, loss, and horizontal gene transfer events. The most ancient event for a given gene represents the first phylogenetic evidence of that gene’s presence, whereas a clustering of events in a particular time interval represents widespread use of the metabolism across diverse lineages. Any lineages which experienced gene events but subsequently became extinct do not leave signatures in the genomic record, so in these cases the origin of each gene could be earlier than we report. Our estimates represent lower bounds.

Estimates made from different clock models vary, with UGAM predicting generally younger events than LN and CIR (Supplementary Figs. 1–14) but all three suggest that microbial communities began importing phosphate before they evolved methods of assimilating reduced phosphorus molecules in the form of phosphonate, phosphite and hypophosphite. For example, the earliest evidence of the phosphate-sodium symporter pnas appears in the Paleoarchean, whilst the earliest evidence for phosphite import and oxidation emerges in the Neoarchean (Fig. 3). For the rest of this paper, we focus on estimates from the CIR clock model because it has been shown to outperform other clock models in terms of model fit43 and consistency between lineages participating in the same HGT events40.

Using this clock model, we conducted further tests to see how ecceTERA’s reconciliation parameters44 impact the estimated histories of phosphorus-cycling enzymes. Duplication, loss, and transfer events involving each phosphorus-cycling gene are predicted based on three things: The shape of the gene tree and molecular clock, the estimated ages of lineages participating in the same transfer events, and the expected cost of transfer in relation to duplication, loss, or speciation of the gene. The latter is a proxy for the general difficulty of duplication in comparison to loss and transfer, so costs remain the same regardless of which genes or species are being assessed. Default costs assume that transfer is more expensive than duplication, loss and speciation in a ratio of 3:2:1:0 because these result in the least genome size variation between parent and daughter lineages45. However, genome size variation is difficult to predict on multi-million-year timescales (Supplementary Discussion), so we conducted reconciliations with different HGT costs (namely 2, 4, and 6) to represent scenarios with more genome size variation45. We found that different costs of horizontal gene transfer (HGT) result in older events for pstS, phnM and phnJ (Supplementary Figs. 15 and 16), but as stated above, microbial communities always begin importing phosphate before assimilating reduced phosphorus molecules, and the earliest evidence for phosphate and phosphite import are in the Paleoarchean and Neoarchean respectively (Supplementary Figs. 15 and 16). For the rest of the manuscript, we report on the findings of all costs tested.

Microorganisms began catabolising phosphonates with the evolution of phnX in the Mesoarchean, but the ability remained confined to a small subset of lineages in comparison to phosphate import mechanisms (Fig. 3, Supplementary Figs. 15 and 16). The gene phnX encodes phosphonatase, which acts on phosphonoacetaldehyde to produce phosphate and small organic molecules for downstream metabolisms within the cell (Fig. 1). As time progressed, the number of potential phosphonate molecules that could be used as a source of phosphorus increased with the emergence of phnW, encoding 2-aminoethylphosphonate-pyruvate transaminase, which can act in concert with phnX to release phosphate from 2-aminoethylphosphonate (2-AEP) (Figs. 3, and 1). This enzyme is estimated to have emerged in either the Mesoarchaean (HGT cost of 6), Neoarchean (HGT cost of 3 or 4) or Paleoproterozoic (HGT cost 2) depending on the event cost (Supplementary Figs. 15 and 16). The genomic repertoire for catabolizing phosphonates also increased in the Paleoproterozoic with the spread of the genes phnJ and phnM, encoding key components of broad-specificity-CP-lyases, which emerged in the Mesoarchean, Neoarchean or start of the Paleoproterozoic depending on event cost (Supplementary Figs. 15 and 16). These genes would have enabled micro-organisms to liberate phosphate and glycine from a post-cursor of 2-AEP and to catabolise a wide variety of different phosphonates molecules (Fig. 3).

Biological methods of producing these phosphonates began to emerge during the Great Oxidation Event (GOE), with the evolution of pepM for transforming the relatively oxidised phosphate ester phosphoenolpyruvate (P in +5 oxidation state) to the more reduced phosphonate phosphonopyruvate (P in +3 oxidation state) and vice versa. This reaction is reversible, so another enzyme is required to force the reaction in the direction of phosphonate production; ppd (otherwise known as aepY) encodes phosphonopyruvate decarboxylase, which is one of the most common methods of doing this19. Here, we find phylogenomic evidence to suggest that ppd emerged less than 300 million years after pepM (Fig. 3, Supplementary Figs. 15 and 16).

The earliest evidence for microbial utilisation of phosphite and hypophosphite appeared in the Neoarchean (or late Mesoarchean if the HGT cost is very high at 6) with the advent of ptxB for importing phosphite and htxA for oxidising hypophosphite and formate (Fig. 3, Supplementary Figs. 15 and 16)42. It has previously been hypothesised that dissimilatory phosphite oxidation (DPO) emerged in the most recent common ancestor of Firmicutes and Deltaproteobacteria25, which emerged ~3.2 Ga46 or ~2.9 Ga47 depending on which molecular clocks are used to infer their age. Here, we predict that this ancestor arose ~3.45 Ga–~3.85 Ga depending on which clock model is used (Supplementary Table 2), but our gene-tree-species-tree-reconciliations do not find evidence of ptxD in this ancestor. Instead, the earliest event indicative of the presence of ptxD is speciation which occurred within the Firmicutes 2.30 Ga (confidence intervals span 2.56–2.22 Ga, Fig. 3). This could range back to 2.83 Ga if HGT was more difficult (Supplementary Fig. 15c), but it places an upper limit on the origin of DPO in the Meso-Neoarchean boundary or more recently unless ptxD was present in an earlier lineage which subsequently went extinct.

Methods
Species tree reconstruction
Our evolutionary tree of life was constructed from information present in the ribosomal proteins of 865 genomes (Supplementary Dataset 1), which were chosen to represent the full diversity of bacteria and archaea using methods described in Mateos, et al.55. To achieve this, we included one representative from each order of the tree of life using information from the Genome Taxonomy Database (GTDB) release 9556,57. Where possible, a representative with the gene ptxD was included to better capture the diversity of phosphite-oxidising organisms. We also include six cyanobacteria, six plastids and four algae to enable us to calibrate the molecular clock using microfossils spanning a range of geological time intervals. From these genomes, a species tree was reconstructed using 16 ribosomal proteins (namely L18, L3, L5, S8, L4, S3_C, L6, L2, L15e, S19, S17, L22, S10, L24, L16 and L14) that have been used previously to infer the tree of life58–60. Different congregations of 14 to 40 proteins34,36,47,58,61–64, sometimes 38165 or 40066, have also been used to reconstruct the tree of life, but most include ribosomal proteins (Supplementary Table 3) because they are ubiquitous and resistant to horizontal gene transfer67. Therefore, homologues of our 16 ribosomal proteins were obtained from each genome with HMMER3 v3.3.268 (95% of the genomes contained homologues of eight or more of these genes), then aligned and trimmed with Muscle v5.169 and TrimAl v1.4.rev1570, all implemented in GToTree v1.6.3471.

The resulting alignments were visualised in AliView v.1.2872 and used to reconstruct a species tree using the maximum likelihood methodology implemented in IQTREE v2.0.373. Partitioned analyses74 were applied to allow each protein to evolve under appropriate substitution models. To find these substitution models, ModelFinder was implemented with default options (no mixture models) and ‘MERGE’75. The best was found to be a single partition containing all proteins, so the LG substitution rate matrix was applied to all 16 ribosomal protein alignments, with 10 categories of the FreeRate model to estimate substitution rates and their variations at different sites. Further variations in substitution rates through time (as described in Lopez et al.76) were modelled by applying the -Q option to account for heterotachy.

Species trees made using different methods sometimes estimate different relationships between phyla (e.g., as demonstrated her64), so we incorporated information from previous studies which thoroughly investigated macroevolutionary relationships between large clades of bacteria and archaea (including 34, 35, 47, 65, 77, 78) by applying a topological constraint on the relationships between Gracilicutes, Terrabacteria, Patescibacteria/CPR and Euryarchaeota (including Halobacteriota, Thermoplasmatota, Methanobacteriota, Hydrothermarchaeota and Hadarchaeota) in our tree (Supplementary Fig. 20). For all analyses, branch supports were estimated with 1000 replicates of ultrafast bootstrap approximations79 and the SH-aLRT test80. Three independent replicates were completed, but only one consensus tree with the best (a.k.a. closest to zero) log-likelihood value was retained for further analyses. The resulting phylogenies were visualised in TreeViewer v2.1.081, FigTree v.1.4.482 and iTOL v6.883.

Divergence time estimation
Our species tree was anchored to a geological timeline by implementing a relaxed Bayesian molecular clock in Phylobayes v.4.184. Topology was fixed to match our species tree, and ages were estimated based on predicted mutation rates for the 16 ribosomal proteins described above (2290 aligned amino acid positions). Substitutions were modelled using a uniform exchangeability matrix and an empirical profile mixture model with 20 profiles (specified using -catfix C20)85. Divergence times were estimated under one of three models to reflect differing views on the inheritance of substitution rates between mother and daughter lineages: These were the autocorrelated lognormal (LN) and Cox-Ingersoll-Ross (CIR) models where substitution rates are inherited from parent lineages43,86, and the uncorrelated gamma multipliers (UGAM) model where substitution rates are chosen independent of ancestry87.

The first diversification in our clock, representing the first radiation of LUCA into bacteria and archaea, was set by applying a uniform distribution spanning 4.4–3.5 Ga to the root node. This is based on the earliest evidence for liquid water on the planet’s surface37,38 and some of the earliest generally accepted microfossils88.

We also implemented eight calibration points based on microfossils of cyanobacteria and algae and geochemical evidence for methanogenesis and oxygenic photosynthesis (Supplementary Table 4). Many of these have been applied previously (e.g., refs. 55, 89–95, but this is the first time they have been combined in a single molecular clock (Supplementary Fig. 21).

For each molecular clock, two independent chains were run for 43,650 to 156,348 cycles, and were considered converged when the relative difference of all parameters was >0.3 and the effective sizes of all parameters were >50. This was assessed using tracecomp implemented in phylobayes v4.184, discarding 10,914 to 49,182 cycles as burn-in. A suitable burn-in was chosen based on traces viewed in Tracer v.1.7.096.

Identification of phosphorus-cycling genes
We searched in our molecular clock genome taxa, for homologs of genes involved in key pathways for the degradation, transport and production of phosphonates, phosphite and hypophosphite using HMMER368. To do this, we used HMM profiles downloaded from the NCBI’s prokaryotic genome annotation pipeline (Supplementary Table 5) and an e-value threshold of 0.1. Where these were unavailable (e.g., for proteins associated with phosphite and hypophosphite transport and oxidation), searches were conducted using the same e-value threshold and custom HMM profiles built from experimentally characterised proteins (Supplementary Table 6) that were aligned in MAFFT v.7.497 with L-INS-I.

The resulting hits were retrieved and aligned with MAFFT v.7.497 using E-large-INS-1. Next, they were trimmed to remove positions comprising >85% gaps with trimAl v1.2 rev5970. This cut-off was chosen because less aggressive trimming has been found to better recover well-established relationships98. Maximum likelihood trees were then generated with IQTREE v.2.0.373 under substitution models chosen by ModelFinder75 with 1000 ultrafast bootstraps79.

To the resulting trees, we annotated bit scores from the HMM searches using the ‘Branchscore’ and ‘Propagate attribute’ modules in TreeViewer v2.0.181. In addition to assigning bit score values to the leaves of the tree, we assigned bit score values to internal nodes based on the average score of their descendants. Any leaves representing homologs which passed the bitscore thresholds of their TIGRFAMs were also highlighted, so that contours could be drawn delimiting clades with the same degree of similarity as the query sequences. Leaves belonging to a monophyletic clade containing the query sequences and any others which passed the respective bitscore threshold were retained for further analyses, whilst the rest were discarded (Supplementary Figs. 22–35). Similar filtering methods have been employed to identify light-harvesting genes91 and arsenic resistance genes62 in prokaryotes.

Phylogenetic analyses of phosphorus-cycling genes
After filtering, the remaining amino acid sequences of each protein family were aligned with MAFFT v. 7.497 and trimmed to remove positions comprising >70% gaps with trimAl70. Bayesian phylogenetic trees were reconstructed using MrBayes v. 3.2.7a99 with a mixed amino acid model prior, invariant sites and gamma-distributed site rates. Chains were considered converged when the average standard deviation of split frequencies was <0.01, the potential scale reduction factor was between 1.00 and 1.02, and the ESS scores of all parameters were >200 after the first 25% of iterations were discarded as burn-in. These were calculated with MrBayes and Tracer v1.796. Bayesian phylogenetic trees of the two phosphate-transporting enzymes, PstS and PNaS, did not converge within three months, so maximum likelihood methodology was implemented instead, using IQ-TREE v2.0.373 with 1,000 ultrafast bootstraps79 and a substitution model chosen by ModelFinder75. To ensure all complex mixture models were included, the following flags were used: -m rate E, I, G, I+G, R and -madd C10, C20, C30, C40, C50, C60, EX2, EX3, EHO, UL2, UL3, EX_EHO, LG4M, LG4X, CF4. Between 100 and 700 iterations were performed before the ultrafast bootstraps converged. A spot check on phosphonatase revealed that similar results were found regardless of whether MrBayes or IQ-TREE were used to reconstruct the evolutionary trees (Supplementary Fig. 36).

Reconciliation of gene trees with the molecular clock
To estimate the timing for duplication, transfer, loss, and speciation events for each gene, each gene tree was reconciled with our time-calibrated species tree using ecceTERA v. 1.2.544. Default settings were applied (event costs: HGT=3, duplication=2, loss=1, speciation=0; transfer to the dead allowed so HGTs can occur with unsampled lineages) with amalgamate=true to find the most parsimonious reconciled gene tree from a set of gene trees representing uncertainty in their topology100. These gene tree sets were obtained from the .t files of the Mr Bayes runs using a custom script which discarded the first 25% as burn-in. From the output of these gene-tree-species-tree reconciliations (specifically the symmetric median reconciliation), we collected the timing of all speciations, transfers, duplications, and losses of each gene by editing Python scripts that were first designed and applied to time the evolution of sulphur-cycling enzymes55. These were applied with Python v. 2.7.5. Each horizontal gene transfer includes two separate events which are counted and output by ecceTERA; a departure of the gene from one lineage and a reception into another lineage. To prevent double-counting and over-inflating the number of transfer events, only reception events are presented.

Phosphate concentration
To investigate whether organisms with low-affinity phosphate transporters and CP-lyases are found in environments with high or low phosphate concentrations, we examined the Ocean Microbial Reference Catalogue v2 (OM-RGC.v2) from the Tara Oceans Project48,49. The OM-RGC.v2 included relative abundances of all clusters of orthologous genes (COG) in 139 Tara Oceans metagenomics samples. We associated each gene of interest with a unique COG accession number using the COG search tool (https://www.ncbi.nlm.nih.gov/research/cog/). To quantify how gene abundance varies with phosphate concentration, regression models of the abundance of each COG were fitted to the raw phosphate concentration using the ‘lm’ function of the ‘stats’ package in R v.4.2.2. All analyses include 129 degrees of freedom.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, the authors present a molecular clock-based analysis of the evolution of phosphorous-metabolizing genes in Bacteria and Archaea. The authors perform phylogenetic and molecular clock analyses of gene and species trees to determine the evolutionary history of these genes and relate their appearance in genomes to predicted amounts of different phosphorous species in the ocean of early Earth. Overall, the topic is of interest to researchers trying to reconstruct the conditions of early Earth and early life’s evolution.\nThe manuscript is well written with an apparently straightforward series of analyses. However, this simplicity comes at a cost as important details of the phylogenetic and molecular clock analyses are missing. My major concerns are:\n1. Phylogenetics: the authors hand-pick ~850 species to reconstruct the phylogeny that they use for all their inferences. While I appreciate the difficulty of estimating large tree of life phylogenies, taxon sampling has been shown to affect the topology of trees, and even more so in the case of prokaryotes where HGT and other factors have such a strong influence on the topology. The cursory overview of the topology estimated by the authors (lines 86-90) does not even scratch the surface of the complexity of prokaryote phylogenies and does not allow to truly evaluate the robustness of this reconstruction. As an example, the authors use the clustering of eukaryotes with Asgard archaea as support to the similarity of their phylogeny to previous ones but this support has a bootstrap of 21, which is very low, even for prokaryotic standards. I was not able to see anywhere, a full list of the species selected, an explanation for their selection, and a detailed phylogeny (either a newick file or a figure that can be zoomed in to look at the details of the whole tree of life).\na. A correlated point is the choice of using 16 ribosomal proteins for the alignment. Also in this case, I empathize with the complexity of the analysis the authors are trying to perform but the use of only a few genes has been shown to affect topology and divergence times. At the very least, a discussion of such potential issues should be present.\n2. Timetree: similarly to the previous point, the authors describe divergence times in a limited way, focusing only on few major time periods (not even relaying specific nodes with their estimated times). The calibrations selected are primarily constrained within the cyanobacteria clade, which is not surprising given they are one of the few groups with information. However, also in this case, there is a large body of literature that discusses the effect of calibration points and their boundaries on estimated times. Yet, there is no mention of this in the manuscript and, apparently, no evaluation of the robustness of their results based on different choices.\n3. Gene/species trees: a fundamental aspect of this analysis is the reconciliation of gene and species trees that has allowed the authors to identify HGT, duplications, and losses and, therefore, estimate the evolution of these genes. I was quite surprised to see that all the information provided effectively mounted to default settings in ecceTERA. The reconciliation of deep gene and species evolutionary history is complex and, surely, affected by the settings chosen in a software, at the very minimum.\nIn summary, I appreciate the effort of the authors and find the topic of great interest. However, I find the description of the study too superficial to truly reflect the complexity of prokaryote evolution. The apparent reliance of the authors on a simplistic methodology for phylogenetic, molecular clock, and gene/species tree reconstruction raises concerns on the validity of the results obtained, as they heavily depend on these three aspects of the analysis.", "reviewer_2": "Boden et al., present a fascinating new study on the history of phosphorus cycling. Their exploration is through molecular clocks where they date the origins of different enzymes and map this history onto chemical and modern oceanographic observations to estimate a qualitative history of phosphorus concentrations.\nThis paper is well organized, written, and presented, however I do have some reservations and suggestions which I detail below.\nFirst, environment is dealt with in a very loose way throughout this manuscript. It is not made clear enough to the reader that local environments could be driving the changes observed without any influence from global seawater. While in parts of the text the authors are careful to note that they cannot identify what environments their high-P or low-P inferences suggest, in other parts of the manuscript they leave the reader with the impression that they are constraining marine phosphate levels. I don’t think that being careful with language here takes away from the manuscript as I find the authors insights important regardless of the local-global nature. Therefore I think stating from the outset the limitations of this approach are important.\nSecond the authors utilize the Tara Oceans data set to explore relationships between their inferences and possible environment. However, the Tara Oceans project only looks at the surface ocean, I think it would be a good idea for the authors to explore other resources available such as the Ocean Protein Portal (https://www.oceanproteinportal.org/) to see if their results display interesting trends at deeper depths.\nThird, I think it is important that the authors explore the role of cell population in their molecular clock analysis. For example, recent works have drawn attention to this point that some uncertainty of molecular clocks may be masked by the assumption of a fixed population size (Crockford et al., 2023). I think this is an important point for the authors to explore.\nFourth, I think the inference regarding BIF deposition is interesting but it is very speculative. Again, I think the paper provides important insights into the history of the P cycle in terms of ordering of events, however, I think attempting to provide constraints on marine phosphorus concentrations or linking findings to the BIF record which itself is mired in controversy, detracts from this work.\nL35 Refs: Laakso and Schrag 2014, 2017, 2019; Derry 2015\nL160 – what is meant by outperform? As the results of such explorations are often untestable with the sedimentary record.\nL228 – but not 1000 times higher than all environments.", "reviewer_3": "In “Timing the Evolution of Phosphorus-Cycling Enzymes Through Geological Time” Boden and co- authors present new estimates on the timing of the origin of specific phosphorus-using enzymes, and then tie these findings to a history of phosphate availability on the early earth. The paper uses phylogenetic analysis to identify the history of microbial genes that use phosphorus in various speciations. They conclude that phosphate has been used for the longest time, followed by phosphonates, then by phosphite and hypophosphite. The findings are novel and of interest to the early-life community, and to those who study biogeochemical element cycling, especially around 2.5 billion years ago.\nAdmittedly I am not wholly qualified to comment on the use of genomic timelines for estimating ages. I’d defer to other reviewers’ expertise to evaluate the methods in the paper. My comments below focus mostly on the conclusions and discussion. I first have some comments for ‘discussion’, which may allow for some expansion of the manuscript, and I highlight a major concern (which may just be a typo) that requires further explanation. I’m happy to recommend publication after some of the discussion and the major concern is addressed.\n***Discussion***\nThe paper presents an interesting finding, which seems to disprove the ‘long half-life’ for phosphite (e.g., Herschy et al. 2018). In other words, the likelihood that phosphite as a major phosphorus species throughout the Archean and into the Proterozoic seems to be disproved. That’s an interesting finding but may be reasonable given that the ‘long half-life’ was based on a solution of phosphite sealed under air that was kept on a bench for 5 years that showed ~0.1% oxidation after that time. Other experiments that argue for a long half-life are based on the slow oxidation of phosphite salts (e.g., CaHPO3) that require high temperatures (~500°C) in order to commence oxidation.\nWith the above noted, shorter half-lives of phosphite may be feasible, for example Gull et al. (Life, 2023) demonstrate that solutions of phosphite when heated to ~80°C in the presence of urea [(NH2)2CO] produce phosphate and polyphosphates, from the oxidation of phosphite. That suggests that there may be other routes to oxidizing phosphite that were not considered in these ‘long half-life’ estimates.\nIn general, the oxidative half-life of hypophosphite appears to be much shorter than phosphite, and therefore the accumulation of hypophosphite in the geologic environment would be uncommon as water is sufficient to oxidize it.\nThe paper is pretty cautious about going further back in time than the mesoArchean (for example, at the origin of life, which makes sense given the tools at hand). Do the authors have any idea on the timeline required to lose a gene if, for example, a component ‘disappears’ from the environment? If phosphite were to diminish to nM or sub-nM concentrations, would the genes to use phosphite still exist if such low-phosphite conditions were present on the early earth for ~100 million years? 500 million years? If such genes would not persist then this may be worth including in the discussion.\nThe near-contemporaneous origin for phosphonate production and phosphite utilization suggests the two may be linked. Indeed, the breakdown of phosphonopyruvate should be sufficient to produce either phosphate or phosphite (both are exergonic reactions). Therefore, once anabolism of phosphonates began in a widespread fashion, phosphite was likely a side-product of this process. This may explain why the two appear somewhat contemporaneous.\n***Major***\nThe description of phosphonate utilization (catabolism and production) at lines 162-182 does not seem to line up with Figure 3. Namely, “Microorganisms began catabolising phosphonates with the evolution of phnX in the Mesoarchean” and “Biological methods of producing these phosphonates began to emerge in the late Neoarchean,”. Figure 3, if the labeling is correct, has pink for production beginning in the MesoArchean, and green for catabolism in the neoArchean. Either the figure caption or the text is in error.\nGiven the textual description, it is likely the colors are mixed up in Figure 3 (or my reading is incorrect). In this case, phosphonate catabolism predates phosphonate anabolism. This still strikes me as unusual. There are sparingly few natural abiotic sources of phosphonates in the environment. The only place they have been found has been in the Murchison meteorite, where they are about 0.1% of the total P content of the meteorite (the remainder is phosphates and phosphides), and an influx of carbonaceous material bearing phosphonic acids seems unlikely. How would catabolism have originated prior to having a source of phosphonates? If production is before catabolism, that would make more sense, but is not the argument presented in the paper.\nSimilar to the above, the paper does link their findings to the estimated abundance of phosphate in the geological environment, and those findings seem consistent. However, for phosphite, phosphonate, and hypophosphite, whether these are presumed to be in the environment or not is less clear. Are the authors arguing that these are trace components in the environment that have always been present, and then when phosphate decreased to such a level (in some environments) that there became an evolutionary pressure to start uptake of these molecules? Or is the manuscript arguing that these molecules began to become more prevalent in the environment at the timing of the evolutionary uptake for these (for example, through biological production)? It seems maybe the former for phosphite and the latter for phosphonate? This could be made clearer. It is easier to envision an abiotic source of phosphite than of phosphonates.\nLine 287 ”In the absence of iron oxidation, all iron would have been in the Fe(II) state, which renders phosphate relatively more soluble, possibly allowing phosphate concentrations of hundreds to thousands of μM”. I think this is probably too strong for the conclusion. For example, Fe2+ with phosphate can form vivianite, which the Brady paper does discuss, but dismisses at the higher pH (~8). Further, one of the most common Fe-minerals is magnetite which is Fe2+(Fe3+)2O4. Magnetite would certainly have been present at this time.\nMinor changes\n---Check the spelling of “phosphorus” at 115.\n---“depleted” at 282\n---Figure 5 is somewhat confusing. The arrows are unclear. Perhaps putting phosphate (genomic estimates) in black but dashed lines may help, as the “geochemical estimates” are also for phosphate."}, "author_response": "Reviewer #1 (Remarks to the Author):\nIn this manuscript, the authors present a molecular clock-based analysis of the evolufion of phosphorous-metabolizing genes in Bacteria and Archaea. The authors perform phylogenefic and molecular clock analyses of gene and species trees to determine the evolufionary history of these genes and relate their appearance in genomes to predicted amounts of different phosphorous species in the ocean of early Earth. Overall, the topic is of interest to researchers trying to reconstruct the condifions of early Earth and early life’s evolufion.\nThe manuscript is well wriften with an apparently straighfforward series of analyses. However, this simplicity comes at a cost as important details of the phylogenefic and molecular clock analyses are missing. My major concerns are:\n1. Phylogenefics: the authors hand-pick ~850 species to reconstruct the phylogeny that they use for all their inferences. While I appreciate the difficulty of esfimafing large tree of life phylogenies, taxon sampling has been shown to affect the topology of trees, and even more so in the case of prokaryotes where HGT and other factors have such a strong influence on the topology. The cursory overview of the topology esfimated by the authors (lines 86-90) does not even scratch the surface of the complexity of prokaryote phylogenies and does not allow to truly evaluate the robustness of this reconstrucfion. As an example, the authors use the clustering of eukaryotes with Asgard archaea as support to the similarity of their phylogeny to previous ones but this support has a bootstrap of 21, which is very low, even for prokaryofic standards. I was not able to see anywhere, a full list of the species selected, an explanafion for their selecfion, and a detailed phylogeny (either a newick file or a figure that can be zoomed in to look at the details of the whole tree of life).\nA full list of the taxa and a newick file of our species tree has been added to our open science framework repository here:\nhftps://osf.io/vt5rw/?view_only=b13a53f4d87c44d1a82a18b176523c5b. They are enfitled ‘TaxaList.csv’ and ‘Q2_merged_constraint3.treefile’.\nWe have also expanded on the relafionships between phyla on lines 87 to 92 to help readers evaluate topology.\nAn explanafion of why we chose the 865 taxa was presented in the methods secfion under ‘Species Tree Reconstrucfion’. We were aftempfing to sample across the diversity of life (as currently understood), but are aware that other methods of sample selecfion exist (e.g. the prototype selecfion algorithm applied in (Zhu et al., 2019) and the quanfile-based equafions of (Marfinez-Gufierrez and Aylward, 2021)). One genome per order with the addifion of 12 cyanobacteria and plasfids for calibrafion purposes were chosen to balance broad and diverse sampling with computafional restricfions.\na. A correlated point is the choice of using 16 ribosomal proteins for the alignment. Also in this case, I empathize with the complexity of the analysis the authors are p0 trying to perform but the use of only a few genes has been shown to affect topology and divergence fimes. At the very least, a discussion of such potenfial issues should be present.\nA table has been compiled summarising the number of proteins used in other tree of life studies (Supplementary Informafion, Table 2 below). Although a couple have used in the region of 300-400 proteins, most use 14-30 to reconstruct the tree of life and esfimate divergence fimes. The excepfion is Zhu et al. (2019), where a dataset of 381 proteins were used to reconstruct the tree of life, but these had to be down- sampled to an alignment of 5,000 amino acid sites (equivalent to ~16 proteins if prokaryofic genes span 924 nucleofides (Xu et al., 2006)) to run the molecular clock.\nSome extra sentences have been added to the methods for ‘Species Tree Reconstrucfion’ to explain that other methods of reconstrucfing phylogenies are available.\nTable 2: The ribosomal proteins used here to reconstruct the tree of life and molecular clock have been utilised in several independent studies for a similar purpose. A complete representation of papers reporting on the topology of the tree of life, including Bacteria and Archaea is included alongside the presence or absence of the ribosomal proteins used in this analyses in their multi-gene datasets. The total number of proteins used in each dataset is listed on the right alongside a column detailing whether they also ran a molecular clock using the same proteins.\n2. Timetree: similarly to the previous point, the authors describe divergence fimes in a limited way, focusing only on few major fime periods (not even relaying specific nodes with their esfimated fimes). The calibrafions selected are primarily constrained within the cyanobacteria clade, which is not surprising given they are one of the few groups with informafion. However, also in this case, there is a large body of literature that discusses the effect of calibrafion points and their boundaries on esfimated fimes. Yet, there is no menfion of this in the manuscript and, apparently, no evaluafion of the robustness of their results based on different choices.\nWe understand your concern. Calibrafion points are key determining factors for any molecular clock study, so it’s important to choose values which are accurate. We thought carefully about which calibrafion points to apply, and chose fossil and geochemical evidence that has been verified and published. For example, our calibrafions on the origin of oxygenic photosynthesis, akinetes and cyanobacterial symbionts have been applied in mulfiple papers by Sanchez-Baracaldo et al (e.g. (Sanchez-Baracaldo, 2015, Sanchez-Baracaldo et al., 2005, Sanchez-Baracaldo et al., 2017, Sanchez-Baracaldo et al., 2014) and our calibrafions for the origin of red algae, eukaryotes and methanogenesis were applied in (Mateos et al., 2023, Parsons et al., 2021). Previous work has tested different age constraints for these lafter three calibrafions based on broadly accepted values in the literature (so-called “conservafive” calibrafion points) and the earliest available evidence (so-called “liberal”)(Parsons et al., p1 2021, Mateos et al., 2023). The liberal points yielded unrealisfic ages for the origin of LUCA, before the origin of planet Earth, so these were discarded and the conservafive points were applied here.\nRecent geochemical evidence on the origin of methanogenesis suggests that it may have evolved earlier than the >2.7 Ga calibrafion we inifially applied. To invesfigate this, we moved this calibrafion earlier to >3.46 Ga to reflect geological evidence of methanogenesis 3.42 Ga from the Barberton Greenstone Belt (Cavalazzi et al., 2021) and 3.46 Ga in the Pilbara Craton (Ueno et al., 2006) as well as phylogenefic evidence derived from horizontal gene transfers which indicate methanogenic Euryarchaeota 3.51 Ga (Wolfe and Fournier, 2018). We ran this with the lognormal model and found that it increases the age of crown Archaea from 3.26 Ga (confidence intervals span 3.45 to 3.08 Ga) to 3.74 Ga (confidence intervals span 3.83 to 3.66 Ga). Few phosphorus-cycling genes are found within Archaea (Figure 2), and the older calibrafion has liftle impact on crown bacteria, so when we reconcile our protein trees with this new clock, there is liftle impact on the fiming of evolufion of phosphorus-cycling enzymes (Figure below).\np2\nFigure: Calibrafing methanogenesis to evolve earlier in Earth history has liftle impact on the esfimated fiming of the evolufion and spread of phosphorus-cycling enzymes. One of seven calibrafion points is altered to force methanogenesis to evolve in the MRCA of TACK and Euryarchaeota > 3.46 Ga instead of > 2.7 Ga to reflect the findings of (Ueno et al., 2006, Wolfe and Fournier, 2018, Cavalazzi et al., 2021). These reconciliafions were calculated using results of molecular clocks made with the lognormal clock model. Mya, million years ago; GOE, great oxygenafion event.\n3. Gene/species trees: a fundamental aspect of this analysis is the reconciliafion of gene and species trees that has allowed the authors to idenfify HGT, duplicafions, and losses and, therefore, esfimate the evolufion of these genes. I was quite surprised to see that all the informafion provided effecfively mounted to default seftings in ecceTERA. The reconciliafion of deep gene and species evolufionary history is complex and, surely, affected by the seftings chosen in a software, at the very minimum.\nWe appreciate that default seftings can reflect a lack of knowledge from the authors, but this was not the case here. We amalgamated our gene trees to account for uncertainty in their topology and roofing and transfer to the dead was allowed to mifigate the impact of reduced taxon sampling (many orders, families, genera and species could have received genes via lateral gene transfer despite not being represented in the fime-calibrated tree of life). While these are well-founded methodological decisions, the relafive cost of transfer vs. duplicafion vs. loss vs. speciafion can be debated, so in response to these concerns we have modified these parameters to assess how they impact our conclusions.\nThe default costs for gene transfers, duplicafions, losses and speciafions have been found to minimise genome size variafion between parent and daughter lineages (David and Alm, 2011). Whether genome size would have remained constant between parent and daughter lineages on fimescales of the order of hundreds of millions of years can be debated. Some findings suggest that genome size distribufions are strongly linked to evolufionary history (Marfinez-Gufierrez and Aylward, 2022) implying liftle change across lineages. In contrast, others have found that terrestrial prokaryotes have larger genomes than aquafic and host-associated prokaryotes, implying that related organisms from different habitats would have different genome sizes (Rodríguez-Gijón et al., 2021).\nIf genome size does vary between parent and daughter lineages that span millions of years, and David and Alm’s findings are valid, then HGT costs would have been higher (David and Alm, 2011). Our results suggest that when HGT is more expensive (cosfing 4 or 6 units instead of 3), phosphorus-cycling enzymes evolve in the same general order with mechanisms of phosphate import emerging in the Paleoarchean, followed by mechanisms of phosphonate catabolism in the Mesoarchean, followed by phosphite ufilisafion and phosphonate producfion in the Neoarchean and during the great oxygenafion event (see figure below). The main differences lie in the fiming of evolufion of the high-affinity phosphate uptake mechanism, pstS. When HGT is more costly, pstS is esfimated to have evolved around the same fime as pnas, not later as default costs imply. Therefore, microbes may have lived in environments with low concentrafions of phosphorus for longer, since the Paleoarchean. Similarly, the fiming of emergence of environments with very low phosphate concentrafions could have been earlier because phnJ and phnM are esfimated to evolve p3 before the GOE, as opposed to after it or towards the end. Linking all scenarios is a lag spanning a few hundred million years between the emergence and spread of pstS (indicafing low phosphate concentrafions) and phnJ and phnM (indicafing very low phosphate concentrafions).\np4\nFigure: Effect of higher HGT costs on the esfimated fiming of the evolufion and spread of phosphorus-cycling enzymes. HGT cost is varied from the default of 3 (a) to 4 (b) and 6 (c) whilst retaining the default costs of gene duplicafion, loss and speciafion. These p5 reconciliafions were calculated using results of molecular clocks made with the CIR clock model. Mya, million years ago; GOE, great oxygenafion event.\nIf HGT was cheaper, cosfing 2 instead of 3 units whilst all other events remain the same, the results remain broadly similar to those with default parameters, although we cannot model the evolufion of pstS because the number of opfimal reconciliafion is too large for ecceTERA to compute (figure below). We have modified the discussion and included these in the supplementary informafion to demonstrate the uncertainty.\nFigure: Effect of lower HGT costs on the esfimated fiming of the evolufion and spread of phosphorus-cycling enzymes. HGT cost is varied from the default of 3 (a) down to 2 (b) whilst retaining the default costs of gene duplicafion, loss and speciafion. These reconciliafions were calculated using results of molecular clocks made with the CIR clock model. Mya, million years ago; GOE, great oxygenafion event.\nThe manuscript has been edited to include the results of these runs.\np6\nIn summary, I appreciate the effort of the authors and find the topic of great interest. However, I find the descripfion of the study too superficial to truly reflect the complexity of prokaryote evolufion. The apparent reliance of the authors on a simplisfic methodology for phylogenefic, molecular clock, and gene/species tree reconstrucfion raises concerns on the validity of the results obtained, as they heavily depend on these three aspects of the analysis.\nReviewer #2 (Remarks to the Author):\nBoden et al., present a fascinafing new study on the history of phosphorus cycling. Their explorafion is through molecular clocks where they date the origins of different enzymes and map this history onto chemical and modern oceanographic observafions to esfimate a qualitafive history of phosphorus concentrafions.\nThis paper is well organized, wriften, and presented, however I do have some reservafions and suggesfions which I detail below.\nFirst, environment is dealt with in a very loose way throughout this manuscript. It is not made clear enough to the reader that local environments could be driving the changes observed without any influence from global seawater. While in parts of the text the authors are careful to note that they cannot idenfify what environments their high-P or low-P inferences suggest, in other parts of the manuscript they leave the reader with the impression that they are constraining marine phosphate levels. I don’t think that being careful with language here takes away from the manuscript as I find the authors insights important regardless of the local-global nature. Therefore I think stafing from the outset the limitafions of this approach are important.\nFair enough, we’ve modified our wording to help clarify this in the discussion.\nSecond the authors ufilize the Tara Oceans data set to explore relafionships between their inferences and possible environment. However, the Tara Oceans project only looks at the surface ocean, I think it would be a good idea for the authors to explore other resources available such as the Ocean Protein Portal (hftps://www.oceanproteinportal.org/) to see if their results display interesfing trends at deeper depths.\nCool resource, thanks! The Ocean Gene Atlas (hftps://tara-oceans.mio.osupytheas.fr/ocean-gene- atlas/) creates some nice graphs of gene abundance with depth, but seems to be limited to a smaller dataset. There are no records of any of the phosphonate cycling enzymes (namely PhnJ, PhnM, PhnZ, PhnX, PhnA, PhnW, Ppd or PepM), or the phosphite cycling enzymes (PtxB and PtxD) or HtxB or PNaS in any of their records from the surface ocean to 6,000m depth, but there are records of PstS, PNaS and HtxA. HtxA is present in a single Marinobacter strain between 50 and 800m depth, whereas PstS is present in 4 strains beneath 200m: Prochlorococcus, Limnobacter, Alcanivorax and Teredinibacter. PNaS is in Desulfococcus up to 1,250m depth. These differences are surprising given the detecfion of some of these genes in the Tara dataset, so may come down to differences in sampling strategies and processing. Previous studies (e.g. (Lockwood et al., 2022)) have looked for correlafions of phosphorus-cycling genes with phosphate concentrafions at depth, and have found liftle correlafion between phosphate concentrafion and the abundance of phosphorus-cycling enzymes at depth. This p7 could be because phosphate concentrafions at depth are higher, so phosphonates are being used as a source of carbon or oxygen instead of phosphorus. Sentences have been added to the discussion to explain that the availability of other molecules could also influence the evolufion of phosphonate catabolism genes.\nThird, I think it is important that the authors explore the role of cell populafion in their molecular clock analysis. For example, recent works have drawn aftenfion to this point that some uncertainty of molecular clocks may be masked by the assumpfion of a fixed populafion size (Crockford et al., 2023). I think this is an important point for the authors to explore.\nCrockford et al (2023) make some good points about how the potenfial for evolufion can change based on the size of the populafion. In a larger populafion, there are more replicafions in a given fime-frame and environment due to the number of cells replicafing. As a result, there is more potenfial for mutafion, so if the number of cells on Earth has increased over geological fime in-line with producfivity (Crockford et al, 2023), the potenfial for mutafion has also increased. This is parfially accounted for in our molecular clocks because they are relaxed, meaning that mutafion rates can vary over fime. These rates vary between lineages, but as yet no method of accounfing for variafion within lineages exist. As a result, the only way to improve resolufion through fime is to include more genomes in the analyses. We were limited to including a few hundred genomes here because adding more genomes to improve resolufion in rate variafion prevents the molecular clocks from converging in a reasonable amount of fime, even on high performance compufing clusters. Therefore, we add a sentence to the discussion to point out this limitafion:\n“The inclusion of novel lineages, or new models of evolufion to account for within-lineage variafion in mutafion rates may alter our conclusions”\nFourth, I think the inference regarding BIF deposifion is interesfing but it is very speculafive. Again, I think the paper provides important insights into the history of the P cycle in terms of ordering of events, however, I think aftempfing to provide constraints on marine phosphorus concentrafions or linking findings to the BIF record which itself is mired in controversy, detracts from this work.\nWe have reworded this point and emphasized that is just speculafive. However, we prefer to keep this suggesfion in the manuscript, because it may offer a potenfial explanafion for the data.\nL35 Refs: Laakso and Schrag 2014, 2017, 2019; Derry 2015\nGreat papers, thanks! Added.\nL160 – what is meant by outperform? As the results of such explorafions are often untestable with the sedimentary record.\nCIR outperforms other clock models (namely LN and UGAM) because it produces befter Bayes factors than UGAM (Lepage et al., 2007) and produces dates which are more consistent than LN and UGAM with relafive dates of different lineages that have parficipated in the same horizontal gene transfers (Fournier et al., 2021). Bayes factors test how well a model of evolufion fits the alignment used to generate a molecular clock, so both of these methods are independent of any geological records. I’ve added some extra words to clarify:\n“we focus on esfimates from the CIR clock model because it has been shown to outperform other clock models in terms of model fit 44 and consistency with relafive dates of lineages parficipafing in the same HGT events 45”\np8\nL228 – but not 1000 fimes higher than all environments.\nAdded the word ‘marine’ to clarify\nReviewer #3 (Remarks to the Author):\nIn “Timing the Evolufion of Phosphorus-Cycling Enzymes Through Geological Time” Boden and co- authors present new esfimates on the fiming of the origin of specific phosphorus-using enzymes, and then fie these findings to a history of phosphate availability on the early earth. The paper uses phylogenefic analysis to idenfify the history of microbial genes that use phosphorus in various speciafions. They conclude that phosphate has been used for the longest fime, followed by phosphonates, then by phosphite and hypophosphite. The findings are novel and of interest to the early-life community, and to those who study biogeochemical element cycling, especially around 2.5 billion years ago.\nAdmiftedly I am not wholly qualified to comment on the use of genomic fimelines for esfimafing ages. I’d defer to other reviewers’ experfise to evaluate the methods in the paper. My comments below focus mostly on the conclusions and discussion. I first have some comments for ‘discussion’, which may allow for some expansion of the manuscript, and I highlight a major concern (which may just be a typo) that requires further explanafion. I’m happy to recommend publicafion after some of the discussion and the major concern is addressed.\n***Discussion***\nThe paper presents an interesfing finding, which seems to disprove the ‘long half-life’ for phosphite (e.g., Herschy et al. 2018). In other words, the likelihood that phosphite as a major phosphorus species throughout the Archean and into the Proterozoic seems to be disproved. That’s an interesfing finding but may be reasonable given that the ‘long half-life’ was based on a solufion of phosphite sealed under air that was kept on a bench for 5 years that showed ~0.1% oxidafion after that fime. Other experiments that argue for a long half-life are based on the slow oxidafion of phosphite salts (e.g., CaHPO3) that require high temperatures (~500°C) in order to commence oxidafion.\nWith the above noted, shorter half-lives of phosphite may be feasible, for example Gull et al. (Life, 2023) demonstrate that solufions of phosphite when heated to ~80°C in the presence of urea [(NH2)2CO] produce phosphate and polyphosphates, from the oxidafion of phosphite. That suggests that there may be other routes to oxidizing phosphite that were not considered in these ‘long half-life’ esfimates.\nIn general, the oxidafive half-life of hypophosphite appears to be much shorter than phosphite, and therefore the accumulafion of hypophosphite in the geologic environment would be uncommon as water is sufficient to oxidize it.\nThe paper is prefty caufious about going further back in fime than the mesoArchean (for example, at the origin of life, which makes sense given the tools at hand). Do the authors have any idea on the fimeline required to lose a gene if, for example, a component ‘disappears’ from the environment? If phosphite were to diminish to nM or sub-nM concentrafions, would the genes to use phosphite sfill exist if such low-phosphite condifions were present on the early earth for ~100 million years? 500 p9 million years? If such genes would not persist then this may be worth including in the discussion.\nHow quickly would a gene disappear if the phosphorus-compound it acts upon disappears from the environment? It’s difficult to know. Genes can appear and disappear from populafions in a few months or years (Bendall et al., 2016, Anderson et al., 2022), but these populafions are relafively small, spanning just one hydrothermal vent or lake. To understand how a gene could be lost from an enfire lineage which probably has representafives in hundreds of populafions in several different geographical areas, would require every place and populafion containing members of that lineage to experience the same. To drive this would require a global loss of the molecules these genes act on. That global loss would also have to persist because genes are often retained in populafions if they are occasionally useful. If the global average dropped from micromolar to sub-nanomolar concentrafions of phosphite, but a few environments (e.g. deep subsurface) retained higher concentrafions and microbial populafions, then the genes for phosphite import and oxidafion would remain. However, if the global average dropped by the same amount, but no phosphite oases existed, the gene would be lost. To predict how quickly a gene would be lost from global biota requires an understanding of habitat heterogeneity, which we can’t predict using the methods in this manuscript.\nThe near-contemporaneous origin for phosphonate producfion and phosphite ufilizafion suggests the two may be linked. Indeed, the breakdown of phosphonopyruvate should be sufficient to produce either phosphate or phosphite (both are exergonic reacfions). Therefore, once anabolism of phosphonates began in a widespread fashion, phosphite was likely a side-product of this process. This may explain why the two appear somewhat contemporaneous.\n***Major***\nThe descripfion of phosphonate ufilizafion (catabolism and producfion) at lines 162-182 does not seem to line up with Figure 3. Namely, “Microorganisms began catabolising phosphonates with the evolufion of phnX in the Mesoarchean” and “Biological methods of producing these phosphonates began to emerge in the late Neoarchean,”. Figure 3, if the labeling is correct, has pink for producfion beginning in the MesoArchean, and green for catabolism in the neoArchean. Either the figure capfion or the text is in error.\nThe figure capfion was incorrect and has now been amended\nGiven the textual descripfion, it is likely the colors are mixed up in Figure 3 (or my reading is incorrect). In this case, phosphonate catabolism predates phosphonate anabolism. This sfill strikes me as unusual. There are sparingly few natural abiofic sources of phosphonates in the environment. The only place they have been found has been in the Murchison meteorite, where they are about 0.1% of the total P content of the meteorite (the remainder is phosphates and phosphides), and an influx of carbonaceous material bearing phosphonic acids seems unlikely. How would catabolism have originated prior to having a source of phosphonates? If producfion is before catabolism, that would make more sense, but is not the argument presented in the paper.\nIt is possible that older genes generafing phosphonates were lost from the genomic record. We have no way of assessing this scenario. Genes that did not survive were likely not widespread across the tree of life at the fime, but they may have existed in high abundance at a local scale, enabling the emergence of phosphonate catabolizing enzymes. Alternafively, it is possible that phosphonates p10 were generated abiofically, perhaps in metamorphic or hydrothermal seftings in the presence of phosphate and organic mafter. To our knowledge, this scenario has so far not been tested experimentally, but it may potenfially offer an explanafion. We have added a brief discussion to this effect to the manuscript:\n“We note that biological phosphonate producfion appears to have emerged after phosphonate ufilizafion in our data, which may indicate an exfinct phosphonate-generafing pathway in the early Archean that was not sufficiently widespread to be preserved in the genomic record. Alternafively, it is possible that minor phosphonates were generated abiofically in metamorphic or hydrothermal seftings in the presence of phosphate and organic mafter, although experimental evidence for this is so far lacking. In any case, our results unambiguously show the Paleoproterozoic radiafion of biological phosphonate uptake pathways, and the mere fact that this signal is preserved in modern organisms speaks to the significance of this metabolic pathway at that fime.”\nSimilar to the above, the paper does link their findings to the esfimated abundance of phosphate in the geological environment, and those findings seem consistent. However, for phosphite, phosphonate, and hypophosphite, whether these are presumed to be in the environment or not is less clear. Are the authors arguing that these are trace components in the environment that have always been present, and then when phosphate decreased to such a level (in some environments) that there became an evolufionary pressure to start uptake of these molecules? Or is the manuscript arguing that these molecules began to become more prevalent in the environment at the fiming of the evolufionary uptake for these (for example, through biological producfion)? It seems maybe the former for phosphite and the lafter for phosphonate? This could be made clearer. It is easier to envision an abiofic source of phosphite than of phosphonates.\nWe can’t differenfiate between these two hypotheses based on genomic data alone. All we can do is look for correlafions between the abundance of a compound and the abundance of a gene to see if the gene is associated with different concentrafions of phosphorus compounds. Some extra sentences have been added to the discussion secfion enfitled ‘Reconstrucfing phosphate availability through fime’ to explain.\nLine 287 ”In the absence of iron oxidafion, all iron would have been in the Fe(II) state, which renders phosphate relafively more soluble, possibly allowing phosphate concentrafions of hundreds to thousands of μM”. I think this is probably too strong for the conclusion. For example, Fe2+ with phosphate can form vivianite, which the Brady paper does discuss, but dismisses at the higher pH (~8). Further, one of the most common Fe-minerals is magnefite which is Fe2+(Fe3+)2O4. Magnefite would certainly have been present at this fime.\nWe softened the wording of this secfion and stress that the BIF explanafion for the declining phosphate trend is somewhat speculafing at this early stage. However, the high solubility proposed by Brady et al. (2022) stands at pH < 8, which is relevant for the Archean (Isson et al., 2018, Krissansen-Tofton et al., 2018, Halevy and Bachan, 2017). Yes, magnefite dates back to 3.8 Ga; however, the total volume of BIF increased in the Neoarchean and Paleoproterozoic, and so it is plausible that oxides were on a global scale a much less important phosphate sink before that fime.\nMinor changes\n---Check the spelling of “phosphorus” at 115.\namended\n---“depleted” at 282\namended\n---Figure 5 is somewhat confusing. The arrows are unclear. Perhaps pufting phosphate (genomic esfimates) in black but dashed lines may help, as the “geochemical esfimates” are also for phosphate.\nwe’ve edited Figure 5 to make it clear that genomic esfimates of phosphate concentrafions are based on different genes."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I would like to thank the authors for the amount of work they added in response to the reviewers' comments. I believe these changes have significantly improved the manuscript and I am happy to recommend it for publication.", "reviewer_2": "I feel the authors have addressed my previous concerns and I am supportive of this MS for publication.", "reviewer_3": "The authors have addressed my points and relabeled the figure that was confusing me. I think the assumptions and findings have been better clarified. I'm happy with this version, but do note that the other reviewers had more comments that are more method-based than mine, and recommend careful consideration of those."}, "author_response": ""}]

Decision letter