Molecular evidence for the early evolution of photosynthetic water oxidation

The evolution of Photosystem II changed the history of life by oxygenating the Earth’s atmosphere. However, there is currently no consensus on when and how oxygenic photosynthesis originated. Here we present a Bayesian relaxed molecular clock analysis of the two core subunits of Photosystem II, D1 and D2, to study the origin of water oxidation. We timed the phylogeny of Type II reaction center proteins using geochemical constraints and calibrations derived from the fossil record of Cyanobacteria and photosynthetic eukaryotes. Age estimates were interpreted in the context of the structure and function of photochemical reaction centers. Firstly, we point out that the ancestral protein to D1 and D2 gave rise to a homodimeric photosystem that was capable of water oxidation, although less efficiently and by a different mechanism than the well-characterized oxygen evolving complex. Secondly, our results indicate that the gene duplication event that led to the divergence of D1 and D2 is likely to have occurred more than a billion years before the last common ancestor of Cyanobacteria. Thirdly, we show that Cyanobacteria did not obtain Type II reaction centers via horizontal gene transfer. Furthermore, the data suggests that the origin of photosynthesis in the early Archaean was necessarily followed by a rapid diversification of all the families of reaction center proteins, which included the divergence of D1 and D2. Our analysis suggests that primordial forms of water oxidation originated relatively soon after the emergence of photochemical reaction centers in the early Archaean.


Introduction
The transition from anoxygenic to oxygenic photosynthesis initiated when an ancestral photochemical reaction center evolved the capacity to oxidize water to oxygen (Rutherford, 1989). Today, water oxidation is catalyzed in the Mn 4 CaO 5 oxygen-evolving cluster of Photosystem II (PSII) of Cyanobacteria and photosynthetic eukaryotes. How and when Type II reaction centers diversified, and how and when one of these reaction centers evolved the capacity to oxidize water are questions that still remain to be answered. While there is agreement that by 3.5 Ga (billion years before the present) a form of anoxygenic photoautotrophy had already evolved (Butterfield, 2015;Nisbet & Fowler, 2014;Tice & Lowe, 2004), the sedimentological and isotopic evidence for the evolution of oxygenic photosynthesis ranges from 3.7 (Frei et al., 2016;Rosing & Frei, 2004) to the Great Oxidation Event (GOE) at ~2.4 Ga . Similarly, molecular clock studies have generated a wider range of age estimates for the origin of Cyanobacteria between 3.5 (Falcon et al., 2010) to 2.0 Ga (Shih et al., 2017). There is thus great uncertainty and no consensus. For this reason, determining when Photosystem II evolved the capacity to oxidize water should greatly advance our understanding of the origin of oxygenic photosynthesis.
The evolution of Type II reaction center proteins has been described and discussed in some detail before (Beanland, 1990;Blankenship, 1992;Cardona, 2015Cardona, , 2016Nitschke & Rutherford, 1991;Rutherford & Nitschke, 1996;Sadekar et al., 2006) and it is presented and schematized in Fig.   1. Type II reaction centers can be divided in two major families: the oxygenic and the anoxygenic Type II reaction centers. The oxygenic Type II reaction center is also known as Photosystem II, and its electron transfer core is made of two homologous reaction center proteins, D1 and D2, exclusively found in Cyanobacteria and photosynthetic eukaryotes. The Mn 4 CaO 5 cluster is located at the Cterminus of D1 (Ferreira et al., 2004). On the other hand, anoxygenic Type II reaction centers are found in phototrophic members of the phyla Proteobacteria, Chloroflexi, and Gemmatimonadetes, with the latter obtaining the reaction center via horizontal gene transfer from a gammaproteobacterium (Zeng et al., 2014). The core subunits of the anoxygenic Type II reaction centers are known as L and M and lack not only an oxygen evolving cluster but any indication that they are involved in the liberation of molecular oxygen.
There is no doubt that D1, D2, L, and M, share a common origin: Beanland (1990) was the first to record this but has been followed by many others (Cardona, 2015;Nitschke & Rutherford, 1991;Sadekar et al., 2006). That is to say that D1, D2, L and M, all descended from a single protein (denoted II in Fig. 1). The earliest event in the evolution of Type II reaction centers can be described as the divergence of this ancestral protein into two new forms, one ancestral to D1 and D2, the oxygenic branch; and a second one ancestral to L and M, the anoxygenic branch (Fig. 1). Hence, D1 and D2 originated from a gene duplication event and together make a monophyletic clade of Type II reaction center proteins, distinct from that which gave rise to L and M (Cardona, 2015(Cardona, , 2016b. The ancestral protein to D1 and D2 will be referred to as D0 and the ancestral protein to L and M will be referred to as K. As a result of the monophyletic relationship of D1 and D2 and the conserved structural and functional characteristics between these two proteins, it is possible to reconstruct traits of the ancestral photosystem. These conserved characteristics, present in both D1 and D2, but absent in L and M, strongly suggest that the ancestral homodimeric photosystem, made of a D0 dimer, had most of the traits known to be associated with the highly oxidizing potential of standard PSII required for water oxidation and associated repair and protection processes (Cardona, 2016;Cardona et al., 2012;Rutherford & Faller, 2003;Rutherford & Nitschke, 1996). Four structural domains differentiate D1 and D2 from L and M ( Fig. 2A). These are: 1. The N-terminus, which is essential for the efficient disassembly and repair of these subunits upon photodamage and occurs by a mechanism that is identical for D1 and D2 Krynicka et al., 2015).

2.
A segment between the 1 st and 2 nd transmembrane helices, which is required for proteinprotein interactions with the small subunits and the extrinsic polypeptides (Cardona, 2015(Cardona, , 2016. The small subunits are required for assembly and stability of the complex and none of them is present in anoxygenic Type II reaction centers (Dobakova et al., 2007;Komenda et al., 2012;Sugiura et al., 2015). The extrinsic polypeptides are important for the stability of the cluster, in particular PsbO (De Las Rivas et al., 2004;Franzen et al., 1985;Roose et al., 2016).
3. An extended loop between the 4 th and the 5 th transmembrane helices required for the coordination of bicarbonate, which is important for regulation and photoprotection (Brinkert et al., 2016;Ferreira et al., 2004).
The strong inference that can be made from these shared structural characteristics is that before the divergence of D1 and D2, the ancestral photosystem was already undergoing evolutionary changes to support water oxidizing chemistry and to deal with singlet oxygen and other reactive oxygen species.
Further evidence for water oxidation before the divergence of D1 and D2 is found at the sequence and cofactor level. A redox tyrosine-histidine pair is conserved in both D1 and D2, this demonstrates that the P/P +• chlorophyll redox pair in this ancestral photosystem was able to generate enough oxidizing power to form a neutral tyrosyl radical on either side of the homodimeric reaction center (Rutherford & Faller, 2003;Rutherford & Nitschke, 1996). In addition, several potential metal ligands are traceable to D0, at least two carboxylate ligands from glutamate residues at positions equivalent to 170 and 189 of standard D1, the carboxylate ligand at the C-terminus, a histidine at position equivalent to 337, and two additional ligands to the Mn 4 CaO 5 cluster provided from the ancestral protein to CP43 and CP47 (Cardona, 2016;Cardona et al., 2015), Fig. 2B. Conserved carotenoids near the peripheral reaction center chlorophylls, Chl Z and Chl D , implies the presence of photoprotective mechanism to quench triplet chlorophyll states and the existence of electron transfer side pathways (Cardona et al., 2012;Krieger-Liszkay et al., 2008), which are in place to prevent the formation of singlet oxygen. Conserved tyrosine residues on both D1 (Y246) and D2 (Y244), which provide hydrogen bonds to the bicarbonate ligand of the non-heme iron (Umena et al., 2011), suggest that this was present in the ancestral homodimeric photosystem, made of D0 (Fig. 2C). The binding of bicarbonate has now been shown to modulate the potential of the quinones as a protective mechanism against singlet oxygen formation (Brinkert et al., 2016). A rational deduction from this is that the ancestral photosystem made of D0 oxidized water accompanied by the unavoidable production of reactive oxygen species, which then required the innovation of protective measures. From that it follows that the heterodimerization process of D1 and D2 was consequently a path to optimization of efficient water oxidation , including photoprotective mechanisms (Krieger-Liszkay et al., 2008), photoactivation of the Mn 4 CaO 5 cluster, assembly, and repair of the protein complex in the presence of oxygen (Nickelsen & Rengstl, 2013;Nixon et al., 2010).
Several types of D1 can be distinguished phylogenetically  and their evolution is schematized in Fig. 1C. The early evolving forms, referred to as atypical D1 forms (G0, G1, G2 in Fig. 1), are characterized by the absence of some of the ligands to the Mn 4 CaO 5 cluster and have been recently found to be involved in the synthesis of chlorophyll f (M. Y. Ho et al., 2016) or the inactivation of PSII (Murray, 2012;Wegener et al., 2015). The late evolving forms, referred to as the standard D1 forms, are characterized by a complete set of ligands to the Mn 4 CaO 5 cluster and are the main D1 used for water oxidation. Among the standard forms there are also several types, which have been roughly subdivided in two groups: the microaerobic forms of D1 (G3) and the dominant form of D1 (G4). The microaerobic forms are thought to be expressed only under low-oxygen conditions. The dominant form, G4, is the main D1 used for water oxidation by all Cyanobacteria and photosynthetic eukaryotes. Most Cyanobacteria carry in their genomes an array of different D1 types, yet every strain has at least one dominant form of D1 (G4). Therefore, all Cyanobacteria descended from a common ancestor that already had evolved efficient oxygenic photosynthesis, had a dominant form of D1, and was able to assemble a standard Photosystem II virtually indistinguishable from that of later evolving strains.
Based on the phylogeny of reaction center proteins several stages in the evolution of oxygenic photosynthesis can be envisaged: the earliest of these stages is the divergence of Type I and Type II reaction center proteins (1, Fig. 1B); this is then followed by the divergence of the anoxygenic family (L/M) and the oxygenic family (D1/D2) of Type II reaction center proteins (2), then by the duplication event that led to the divergence of D1 and D2 (3), and the subsequent (7) gene duplication events and specializations that created the known diversity of D1 forms, which ultimately resulted in the emergence of the standard form of D1. If D0 was able to support water oxidation chemistry, then the earliest stages in the evolution of oxygenic photosynthesis occurred between stages 2 and 3 as depicted in Fig. 1B. Therefore, the span of time between D0 and the ancestral standard form of D1 (marked 8 in Fig. 1C) represents the duration of the evolutionary trajectory of Photosystem II from a simpler homodimeric water-splitting reaction center to the more complex enzyme inherited by all organisms capable of oxygenic photosynthesis. We denote this span of time by such as a few million years or less for example, then the evolution of oxygenic photosynthesis was likely a sudden and fast process occurring shortly before the GOE as suggested by some recent analyses (Shih et al., 2017;Ward et al., 2016). On the other hand, if Δ T is large: several hundred million years or more for example, then water oxidation could significantly predate the GOE as suggested by some geochemical Satkoski et al., 2015) and phylogenetic data (Blank & Sanchez-Baracaldo, 2010;Schirrmeister et al., 2013).
Here we report an in-depth evolutionary analysis of Type II reaction center proteins including Bayesian relaxed molecular clocks under various scenarios for the origin of photosynthesis. The data presented here indicate that the origin of light-driven water oxidation started early in the evolutionary history of photosynthesis, likely during the early Archaean and long before the last common ancestor of Cyanobacteria.

Change in sequence identity as a function of time
A first approximation to the evolution of Type II reaction centers as a function of time can be derived from the level of sequence identity between D1 and D2 of different species with known divergence times as shown in Fig. 3. For example, the D1 protein of the dicotyledon Arabidopsis thaliana shares 99.7% amino acid sequence identity with that of the dicotyledon Populus trichocarpa, and these are estimated to have diverged between 127.2 and 82.8 Ma (Clarke et al., 2011), see   (Butterfield, 2000) and recently described fossils could push this date back to 1.6 Ga (Bengtson et al., 2017). At the other end of this evolutionary line, the three dominant forms of D1 from Gloeobacter violaceous share on average 79.2% sequence identity with that of C. merolae or 78.5% with that of A. thaliana. If the percentage of sequence identity between pairs of species is plotted as a function of their divergence time, a linear decrease of identity is observed among reaction center proteins at a rate of less than 1% per 100 million years (supplementary Table S2). The trend in Fig. 3 indicates that the rate of evolution of D1 and D2 since the GOE and since the emergence of photosynthetic eukaryotes has remained relatively slow until the present time, if considered over a large geological time scale, with less than 20% change in sequence identity in the past 2. Let us reiterate that all the evidence suggests that all reaction center proteins originated from a single ancestral protein that diversified as the multiple groups of photosynthetic bacteria arose. As a result of this common ancestry, any standard D1 shares on average about 29% sequence identity with any D2 across their entire sequence. Any standard D1 or D2 share on average 17% sequence identity with any L or M. The level of sequence identity falls well below 5% if any Type II reaction center protein is compared with any Type I reaction center protein (Cardona, 2015). In consequence, the rate of evolution of D1 and D2 since the GOE, as estimated from the decrease of sequence identity (<1% per 0.1 Ga), is too slow to account for the evolution of photochemical reaction centers within a reasonable amount of time (Fig. 3, dashed line). In other words, the rate of evolution of reaction center proteins since the origin of life could not have been constant, and any scenario for the origin of photochemical reaction centers at any point in the Archaean requires initially faster rates of evolution than any rate observed since the Proterozoic (Fig. 3, light blue line).
Taking into consideration that D1 and D2 share only about 29% sequence identity, this would suggest, as illustrated in Fig. 3, that the event that led to the divergence of D1 and D2 is more likely to have occurred closer to the origin of the primordial reaction center proteins during the emergence of the first forms of photoautotrophy in the early Archaean, than closer to the Great Oxidation Event.

Bayesian relaxed molecular clock analysis
While the approach used above indicates that the divergence of D1 and D2 is likely placed well before the GOE, it is possible to make further inferences on the evolution of oxygenic photosynthesis by applying a Bayesian molecular clock approach to the phylogeny of Type II reaction center proteins and from the resulting rates of evolution. Figure 4 shows a Bayesian relaxed log-normal autocorrelated molecular clock built using the CAT + GTR + Γ model allowing for flexible boundaries on the calibration points (Lartillot et al., 2009). As an informed starting point we first specified the age of the root (root prior) at 3.5 Ga with a standard deviation of 0.05 Ga. That is to say, that the most ancestral form of a Type II reaction center protein is assumed to have already evolved by 3.5 Ga. Under these conditions, the last common ancestral protein to the standard form of D1 prior to the divergence of the G3 and G4 types ( To test the effect of different root priors on our results we varied the age of the root and the standard deviation over a broad range. Table 1 lists estimates of divergence times of key ancestral Type II reaction center proteins and the respective Δ T value using different root priors. For example, under the assumption that Type II reaction centers had already evolved by 3.8 Ga ago (Czaja et al., 2013;Nisbet & Fowler, 2014;Rosing, 1999), Δ T is found to be 1.17 Ga. Similarly, if it is assumed to be a late event centered at 3.2 Ga, though unlikely, Δ T is still 0.90 Ga. Furthermore, increasing the standard deviation on the root prior pushes the timing of the earliest events in the evolution of the Type II reaction center to even older ages rather than younger ages. Supplementary Table S3 lists the changes in Δ T by varying the standard deviation in the root prior. For example, a root prior of 3.5 Ga with a standard deviation of 0.1 Ga pushes the estimated time for the root to 3.65 Ga; and with a standard deviation of 0.2 Ga to 4.00 Ga. In the former case, Δ T shifts to 1.16 Ga and in the latter to 1.31 Ga.
In Fig. 4 it can also be seen that the divergence of the L and M subunits occurs after the divergence of D1 and D2. The estimated time for the divergence of L and M is centered at 2.87 ± 0.16 Ga, while the time for the divergence of D1 and D2, as we saw above, is 3.22 ± 0.19 Ga. The late divergence of L and M suggests that the divergence of D1 and D2 predates significantly the diversification event that led to the evolution of Chloroflexi-type and Proteobacteria-type L and M subunits. These results highlight the conclusion that it is unlikely that Cyanobacteria obtained Type II reaction centers via horizontal gene transfer (HGT) of L and M from an early evolving representative of the phylum Chloroflexi or Proteobacteria (Cardona, 2015).
The Bayesian clock using flexible boundaries on the calibration points consistently produced ages for the divergence of the D2 subunit of G. violaceous and the dominant form of D1 (G4) after the GOE, in agreement with the work by Shih et al. (2017). Yet, previous molecular clocks have suggested that the last common ancestor of Cyanobacteria might predate the GOE (Schirrmeister et al., 2015), so we also performed a similar analysis that allowed us to explore this scenario. This was achieved using an empirical amino acid substitution model (LG + Γ ) instead of the non-parametric approach described above. We found this to be the only way to locate the D2 of G. violaceous and the dominant form of D1 (G4) before the GOE. The effect of less flexible boundaries on the estimated divergence times is shown in Table 1. Assuming a root prior of 3.5 ± 0.05 Ga, the estimated divergence time for the standard form of D1 becomes 2.64 ± 0.15 Ga and pushes D0 back to 3.40 ± 0.09 Ga, making Δ T 0.77 Ga. On the other hand, if we allowed flexibility on the root prior by increasing the standard deviation to 0.4 Ga (supplementary Table S3), the estimated divergence time for the standard form of D1 becomes 2.83 ± 0.21, but the estimated age of the root is pushed back to 3.92 ± 0.21 Ga with D0 at 3.79 ± 0.21, making Δ T 0.96 Ga. Overall, placing the last common ancestor of Cyanobacteria before the GOE pushes the gene duplication event that led to the divergence of D1 and D2 even closer to the origin of Type II reaction centers and to the origin of photosynthesis.

Rates of evolution
The inferences derived from Fig. 3 reveal that the rates of evolution had to be faster in the initial stages during the Archaean compared with the Proterozoic, even when Δ T is as large as 1.0 Ga. To gain a better understanding of the changes of the rate of evolution of Type II reaction center proteins we plotted the rates as a function of divergence time. In Fig. 5A the rate of evolution (ν) of each node in the tree, expressed as amino acid substitutions per site per unit of time, is plotted against the estimated divergence time for each respective node. It can be seen that the rate at the earliest stage is much faster than the rates observed since the Proterozoic. Clearly, faster rates are necessary to explain the origin and evolution of Type II reaction centers at any point in the Archaean and regardless of when exactly photosynthesis originated, as seen in Fig. 5B. The decrease in the rate of evolution is consistent with the observations derived from Fig. 3 and can be roughly fitted with a first-order exponential decay curve (fitting parameters are presented in supplementary Table S4). Fig. 5A additionally shows that L and M have been evolving at a slightly faster rate than D1 and D2. From this slow-down of the rates of evolution it can be calculated that since each respective duplication event (stages 2 and 4 in Fig. 1) it took about 168 million years for D1 and D2 to fall to 50% sequence identity and about 115 million years for the same to occur for L and M.
The maximum rate of evolution in the D1 and D2 family of reaction center proteins is placed at the node that represents D0. We will refer to this rate as ν max . Fig. 5A shows that the rate of evolution flattens out to comparatively slow rates during the Proterozoic. These rates correspond to the rates of Group 4 D1 and D2. We will refer to the average rate of evolution during this zone of slow change as ν min and it is calculated as the average rate from each node in Group 4 D1 and D2. In substitutions per site per Ga. Based on these assumptions it is almost certain that reaction centers had evolved by 3.5 Ga, to account for the divergence of D1 and D2 in one billion years (ΔT), since the initial rate of evolution had to be about 40 times larger than that observed since the last common ancestor of Cyanobacteria. Table 2 lists the rates of evolution of a diverse number of proteins reported in independent studies in a broad range of organisms. It is found that the rate of evolution of the core subunits of PSII is similar to the rate of other proteins that are very ancient and highly conserved such as membrane proteins involved in bioenergetics or the ribosome. Our estimated rates fall well within the expected range of other cyanobacterial proteins, thus validating our calibration choices and consistent with the expected level of sequence identity as derived from Fig. 3.
Our Bayesian analyses suggests that the evolution and diversification of PSII is better explained by a long span of time since the gene duplication even of the homodimeric core until the emergence of standard PSII inherited by crown group Cyanobacteria. Notwithstanding, a fast rate of evolution at the earliest stage implies that ν max would increase if Δ T is considered to be smaller, as would be the case for an evolutionary scenario in which oxygenic photosynthesis evolves rapidly before the GOE after events of horizontal gene transfer. We illustrate this concept in Fig. 5B  per site per Ga, which is greater than the rate of evolution of some extant viruses, see Table 2.

Sensitivity analysis
To test the reliability of the method we explored a range of contrasting models. We compared the effect of the model of relative exchange rates on the estimated divergence times: Supplementary Fig.   S1 provides a comparison of estimated divergence times calculated with the CAT + GTR model (Z. Yang & Rannala, 2006) against divergence times calculated using the CAT model with a uniform (Poisson) model of equilibrium frequencies. The GTR model does not have a strong effect on the calculated divergence times as the slope of the graph does not deviate from unity when paired with the uniform model (see supplementary Table S6 for linear regressions). Thus under a root prior of 3.5 ± 0.05 Ga a CAT + Poisson model also generated a Δ T of 1.02 Ga, see Table 4.
To understand the effects of the oldest calibration point (point 11, Fig. 4) on the estimated divergence time, we tested a second set of boundaries restricting this point to a minimum age of 2.7 Ga (Calibration 2) under the assumption that the divergence of the genus Gloeobacter may significantly predate the GOE. Supplementary Fig. S2 provides a comparison of the two calibration choices on the overall estimated times for flexible and non-flexible models. If the divergence times using both calibrations are plotted against each other, a linear relationship is obtained (supplementary Table S7). Calibration 2 did not seem to have a very strong effect on the estimated divergence times nor Δ T. For example, under a root prior of 3.5 ± 0.05 Ga and employing Calibration 2, Δ T was 0.97 Ga, see Table 4. We also tested the effect of removing the oldest calibration point from the analysis (point 11, Fig. 4), but this only shifted Δ T to values greater than 1.0 Ga.
In contrast, the choice of model for the evolution of substitution rates had a strong impact on the estimated divergence times as shown in supplementary Fig. S3 and Table S8. Supplementary Fig.   S3 presents a comparison of divergence time estimates of a tree calculated using a relaxed log-normal autocorrelated molecular clock with a tree calculated using an uncorrelated gamma model on the rates of evolution. The autocorrelated model assumes that neighboring branches are more likely to evolve at a similar rate, while the uncorrelated model assumes that the rate of evolution of each branch can vary independently (S. Y. W. Ho & Duchene, 2014;Lepage et al., 2007). Under the uncorrelated model the estimated divergence times of many nodes were aberrantly shifted to younger ages: for example, most Cyanobacterial and eukaryotic D1 clustered in the range of 700 to 0 Ma, which is inconsistent with the fossil record. The molecular mechanism behind this difference could be related to the fact that photochemistry imposes a strong constraint on the evolution of reaction center proteins: as all of them must coordinate and maintain all redox cofactors, chlorophylls, quinones, carotenoids, and a non-heme iron, at a precise orientation and distance from each other to allow for control of electron transfer rates and redox potentials. These rates and potentials are crucial not only for function but also for protection against the formation of reactive oxygen species (Cardona et al., 2012;Rutherford et al., 2012). It seems reasonable then, that the rates of evolution of all Type II reaction center proteins should be correlated between closely related groups, thus corresponding to an autocorrelated model.

Change in sequence identity as a function of time
As an approximation to the evolution of the core subunits of PSII we plotted the level of sequence identity of D1 and D2 as a function of known divergence times (Fig. 3). Two main conclusions can be derived from this plot, which are also in agreement with the molecular clock analysis. Firstly, the three earliest stages in the evolution of photochemical reaction centers: the divergence of Type I and Type II reaction centers, the divergence of anoxygenic and oxygenic families of Type II reaction center proteins, and the divergence of D1 and D2, are more likely to have occurred soon after the origin of the first reaction centers rather than near the GOE. Taking into consideration that there is strong evidence for photosynthesis at 3.5 Ga (Butterfield, 2015;Nisbet & Fowler, 2014;Tice & Lowe, 2004), then all three stages could well predate this time. Secondly, the last common ancestor of crown group Cyanobacteria is more likely to have lived near the time of the GOE rather than shortly after the origin of photosynthesis in the early Archaean. This common ancestor must have had a standard PSII, which places it relatively far after the origin of photosynthesis. The early divergence of D1 and D2 means that the earliest stages in the evolution of photosynthetic water oxidation predates the last common ancestor of Cyanobacteria by about 1.0 Ga.

Bayesian relaxed molecular clock analysis
The application of a Bayesian molecular clock analysis to the phylogeny of Type II reaction centers is strongly constrained by three pieces of well-supported evidence: 1) evidence of photosynthesis at 3.5 Ga, 2) by strong evidence for the GOE at 2.4 Ga, and 3) by the relatively slow rate of evolution of the core proteins over the last 2.0 Ga. Under these constraints, the divergence between D1 and D2 is better explained by water oxidation starting early in the evolutionary history of photosynthesis, in the early Archaean, with the appearance of standard PSII occurring after a long period of evolutionary innovation. We highlight this long period by introducing the concept of Δ T, which defines two specific points in time in the evolution of PSII: the time from the gene duplication event of the homodimeric core to the ancestral standard heterodimeric PSII retained by extant Cyanobacteria. The magnitude of Δ T is dictated by the large phylogenetic distance between D1 and D2 and the relatively slow rate of evolution determined from the geochemical and fossil record. Therefore, it is not surprising that under most models employed in this analysis Δ T is always in the range of 0.7 to 1.2 Ga.
We have considered two possible evolutionary scenarios that are both consistent with a large Δ T (Fig. 6). In the first scenario the standard forms of D1 start to diverge at about 2.4 Ga, as seen in Fig. 4, and diversify into G3 and G4 after the GOE. If we consider that the last common ancestor of Cyanobacteria had a G4 D1, this would set it after the GOE. This scenario, derived from the application of a relaxed molecular clock using a non-parametric CAT model with flexible boundaries, is in agreement with the recent observations by Shih et al. (2017) and other molecular clock studies that placed the divergence of Gloeobacter after the GOE (David & Alm, 2011;Feng et al., 1997;Marin et al., 2017). In this scenario, assuming that the earliest events in the history of photosynthesis started about 3.5 Ga, the divergence of D1 and D2 is set at about 3.2 Ga.
In the second scenario, we considered that the last common ancestor of Cyanobacteria occurs before the GOE as suggested by other molecular clock analyses (Falcon et al., 2010;Sanchez-Baracaldo, 2015;Schirrmeister et al., 2015). In the present work, this scenario can be fitted most satisfactorily with the application of a relaxed molecular clock using an empirical amino acid substitution model (LG + Γ ). In this scenario, under a root prior of 3.5 Ga, the appearance of the ancestral standard form of D1 is set at about 2.6 Ga; and this has the consequence of pushing the divergence of D1 and D2 closer to the root, and thus D0 is set at about 3.4 Ga (Table 1, Fig. 6). This effect is due in part to the fact that the phylogenetic distance between D1 and D2 is invariable, and thus under any scenario the data is better explained by a long span of time separating D0 and the standard heterodimeric PSII. What can be concluded from this is that the older the last common ancestor of Cyanobacteria is, the more likely it is that the divergence of D1 and D2 occurred near the origin of photochemical reaction centers and thus, near the origin of photosynthesis.

Recently, Shih et al. (2017) reported a phylogenetic analysis of Cyanobacteria and their
closest non-photosynthetic relatives, the Melainabacteria (Soo et al., 2017;Soo et al., 2014). Shih et al. (2017) postulated that Cyanobacteria acquired photosynthetic capacity via horizontal gene transfer after they diverged from the Melainabacteria. They hypothesized that the transfer of photosynthetic genes occurred from anoxygenic photosynthetic bacteria into a non-photosynthetic Cyanobacterial ancestor; with water oxidation evolving quickly before the GOE. This scenario requires that the transferred Type II reaction center was still homodimeric occurring before the divergence of D1 and D2. The phylogeny of Type II reaction centers and our Bayesian analysis do not lend support to this scenario because the phylogenetic distance between D1 and D2 is too great and the rates of evolution of the core subunits are too slow for this transfer to have occurred not long before the GOE. Our analysis strengthens the notion that early-branching non-photosynthetic Cyanobacteria emerged after losing water oxidation catalysis in a process of adaptation to heterotrophic lifestyles. The loss of oxygenic photosynthesis is a relatively common process, several well-known examples are: the secondary loss of all PSII genes by the symbiont Atelocyanobacterium talassa (Cornejo-Castillo et al., 2016;Thompson et al., 2012), the complete loss of photosynthesis of the cyanobacerium endosymbiont of Epithemia turgida (Nakayama et al., 2014), the case of parasitic Apicomplexa (Moore et al., 2008) and holoparasitic angiosperms (Ravin et al., 2016).

Rates of evolution
It should be observed that if a relatively constant rate of evolution were to be applied to the phylogeny of Type II reaction centers, the root of the tree would be placed long before the formation of the planet (dashed red line in Fig. 3). Given the fact that the rate of evolution of D1 and D2 is constrained at comparatively slow rates from the Proterozoic until present times, the period of fast evolutionary change has to be located at the earliest stages to account for the evolution of photosynthesis within a reasonable timeframe. We found that even considering an origin of Type II reaction centers at 3.8 Ga the maximum rate during the evolution of D1 and D2, ν max , is many times greater than the measured rates in the Proterozoic.
We are not aware of any systematic study of the change in the rate of evolution of proteins throughout geological time and such a study is outside the scope of this article. However, there are several families of proteins that are likely to follow a similar trend. For example, the phylogeny of rubisco is characterized by limited sequence divergence at the phylum level, but a much greater phylogenetic distance between orthologous groups (Kacar et al., 2017). Indeed, Kacar et al. (2017) suggested that the divergence of Group 3 rubisco of methanogenic Archaea and bacterial photosynthetic Group 1 rubisco occurred as the latter acquired adaptations to deal with oxygen and more oxidizing conditions. The divergence of Group 3 and Group 1 Rubisco also requires initially higher rates of evolution. A similar situation is expected for methanogenic and photosynthetic phosphoribulokinases, which also show a remarkable phylogenetic distance (Kono et al., 2017).
Another similar case is found between the homologous enzymes Ni-sirohydrochlorin a,c-diamide reductive cyclase of methanogeneic Archaea, protochlorophyllide and clorophyllide reductase of phototrophic bacteria (Hu & Ribbe, 2015;Zheng et al., 2016), and nitrogenases. These enzymes show a very high level of sequence identity within phyla due to slow rates of evolution, but a great phylogenetic distance between orthologous groups.
The phenomenon of an initial fast rate of evolution followed by an exponential decrease demands an explanation. We can think of two mechanisms that may account for this observation, but there may be more. The first one is the temperature dependent deamination of cytosine, as suggested by Lewis and coworkers (2016). They calculated that as the Earth cooled during its first 4 Ga, the rate of spontaneous mutation would have fallen exponentially by a factor of more than 4000. That is to say that the rate of spontaneous mutation during the earliest stages in the history of life would have been about three orders of magnitude greater than those observed since the Proterozoic. Lewis and coworkers (2016) calculated that 50% of all spontaneous mutations occurred in the first 0.2 Ga, which matches well with the exponential decay trend seen in Fig. 5A, especially if an origin of photosynthesis is considered to be about 3.8 Ga. The second possibility is higher UV radiation on the planet's surface during the early Earth in the absence of an ozone layer, which could have resulted in rates of DNA damage up to three orders of magnitude greater than in present day Earth, as calculated by Cockell (2000); this higher rate of damage may have led as well to faster rates of change.
Alternatively, both mechanisms could have contributed simultaneously.

Heterodimerization of ancestral Type II reaction centers
Another finding of the molecular clock analysis is that the divergence of L and M seems to occur after the divergence of D1 and D2. This is due to the slightly faster rates of evolution observed for L and M at all periods of time. If we add to this the fact that L and M subunits are about 50 to 90 amino acids shorter than D1 and D2, it should take less time for L and M to converge into the ancestral node K.
This raises the following question: if the heterodimerization of D1 and D2 was driven by optimization of water oxidation and the incorporation of photoprotective mechanisms, what selective pressures triggered the heterodimerization of L and M?
It has been proposed that the heterodimerization of Type II reaction centers was driven by the optimization of electron transfer efficiency and the avoidance of photodamage (Cardona et al., 2012;Rutherford & Faller, 2003;Rutherford & Nitschke, 1996;Rutherford et al., 2012). In a homodimeric Type II ancestor containing two identical exchangeable quinones, charge separation would result in enhanced back-reactions caused by electron transfer to a quinone site when empty or when occupied by a reduced form of the quinone (Rutherford & Faller, 2003;Rutherford & Nitschke, 1996); and additionally, by shorter back-reaction pathways (Cardona et al., 2012;Rutherford et al., 2012). Backreactions from the early radical pairs give rise to chlorophyll triplet states in the heart of the reaction center (Dutton et al., 1972;Rutherford et al., 1981) and these triplet states are dangerous when they encounter oxygen, giving rise to the highly damaging singlet oxygen (Krieger-Liszkay et al., 2008;Rutherford et al., 2012;Vass & Cser, 2009). In the oxygenic D0 homodimeric ancestor of PSII, avoidance of photodamage could be a significant selective pressure for heterodimerization . It seems clear that the enhanced back-reactions intrinsic to the homodimeric Type II reaction centers would have become a major problem whenever oxygen was present. The inefficient D0 reaction center, as the primordial oxygenic photosystem, would have encountered this problem first and thus come under strong selection pressure toward heterodimerization at an early time. The results of the present work fit with this picture. Furthermore, our present finding that heterodimerization started significantly later in the anoxygenic homodimeric Type II reaction center (K) compared to its oxygenic counterpart (D0), can be explained in this context as well, since K would have encountered these selection pressures much later, at a time when oxygen concentrations began to rise in localized environments in the late Archaean period, nearing the GOE (Bosak et al., 2009;Lyons et al., 2014;Riding et al., 2014).
It is rather remarkable that anoxygenic Type II reaction centers contain an asymmetrically located carotenoid in contact with a core bacteriochlorophyll (Deisenhofer & Michel, 1989). The role of this carotenoid is to quench bacteriochlorophyll triplet states in order to prevent the formation of singlet oxygen (Cogdell et al., 2000). In addition, light-harvesting complexes in most anoxygenic photosynthetic bacteria contain carotenoids, which perform a similar photoprotective function (Kim et al., 2007;Melo et al., 2000;Tsukatani et al., 2012). As an extension of this, it does not seem unreasonable to think that even ancestral populations of anoxygenic phototrophic bacteria were under strong selective pressure by the threat of bacteriochlorophyll and chlorophyll triplet-induced formation of radical oxygen species early during the evolutionary history of photosynthesis.

Concluding remarks
The evolution of Type II reaction centers highlights the long history of photosynthetic water oxidation. We showed that the span of time between the gene duplication event that led to D1 and D2 and the appearance of standard PSII could have been about a billion years, placing water oxidation in the early Archaean. So what happened during this long period of time? To begin with, water oxidation originated in a simpler, single-subunit, homodimeric photosystem in a completely anaerobic world. Therefore, the large increase in the structural complexity of PSII, PSI, and associated light harvesting complexes had to occur along this trajectory. This includes the acquisition of many peripheral and extrinsic protein subunits and the heterodimerization of D1 and D2, CP43 and CP47, and PsaA and PsaB. At the same time, the thermodynamic coupling between both photosystems and the retuning of the entire electron transport chain and all electron carriers to increasingly oxidizing conditions also had to occur. This expansion in complexity had to be coupled with the evolution of highly organized assembly and repair processes. Thus, the first water-oxidizing reaction centers may have been active only for briefs amounts of time in the absence of efficient repair. Greater water oxidation efficiency also needed the innovation of photoprotective mechanisms acting at different time scales spanning several orders of magnitude, like dissipatory recombination pathways, non-photochemical quenching, state-transitions, or chromatic adaptation. Furthermore, the light reactions of photosynthesis had to be linked to carbon fixation and other downstream metabolic process. Signaling, feedback, and regulatory mechanisms had to be put in place to control photosynthesis under varying environmental conditions and across a diel cycle. Needless to say, all anaerobic reactions and processes inhibited by oxygen originally found in the earliest anaerobic water-oxidizing ancestors had to be separated from oxygen production or readapted to work under aerobic conditions. The link to carbon fixation is of particular importance. CO 2 levels in the atmosphere were higher than now (Nutman et al., 2017;Sheldon, 2006). However, limiting diffusion across boundary layers in mats and stromatolites would have restricted anoxygenic and early oxygenic phototrophs alike. The development of water oxidation would have opened up the way to faster photosynthetic rates, spurring on gross primary production rates later in the Archean, with the concomitant need for increases in nitrogen fixation.

Phylogeny of Type II reaction centers
Sequences were retrieved from the RefSeq NCBI database using PSI-BLAST restricted to Cyanobacteria, Proteobacteria, and Gemmatimonadetes. A total of 1703 complete sequences were downloaded and aligned using Clustal Omega employing ten combined guide trees and Hidden Markov Model iterations (Sievers et al., 2011). To confirm that the alignment conformed with known structures, the 3D structures of the D1, D2, L and M, from the crystal structures 3WU2 (Umena et al., 2011) of PSII and 2PRC (Lancaster & Michel, 1997) of the anoxyenic Type II reaction center were overlapped using the CEalign (Jia et al., 2004) plug-in for PyMOL (Molecular Graphics System, Version 1.5.0.4 Schrödinger, LLC) and structural homologous positions were cross-checked with the alignment. Maximum Likelihood (ML) phylogenetic analysis was performed using PhyML 3.1 (Guindon et al., 2010) using the LG model of amino acid substitution. The amino acid equilibrium frequencies, proportion of invariable sites, and across site rate variation were allowed to be determined by the software from the dataset. Tree search operations were set as the best from the Nearest Neighbor Interchange and Subtree Pruning and Regrafting method, with the tree topology optimized using five random starts. The ML tree using all sequences is shown in Fig. 1 and it replicates earlier evolutionary studies of Type II reaction centers that used simpler methods and fewer sequences (Beanland, 1990;Cardona, 2015).

Change in sequence identity as a function of time
To get a better understanding of the evolutionary trends of D1 and D2 reaction center proteins as a function of time, we compared the percentage of sequence identity of D1 and D2 proteins from species of photosynthetic eukaryotes with known or approximate divergence times. In total 23 pairs of sequences were compared and are listed in Supplementary Table 1. When two sequences were of different length, the longest was taken as 100%. Of these 23 pairs, the first 16 were based on the fossil calibrations recommended by Clarke et al. (2011) after their extensive review of the plant fossil record. Divergence times were taken as the average of the hard minimum and soft maximum fossil ages suggested by the authors. The last 7 comparisons are approximate dates taken from a recent molecular clock analysis of the evolution of red algae (E. C. Yang et al., 2016). The plot of sequence identity vs approximate divergence time was then fitted with a linear regression and the fitting parameters are shown in supplementary Table S2.

Bayesian relaxed molecular clock and fossil calibrations
A total of 54 bacterial sequences spanning the entire diversity of Type II reaction centers were selected for Bayesian molecular clock analysis, including atypical and standard forms of D1, Alpha-, Beta-, Gammaproteobacteria, Chloroflexales, and Gemmatimonas phototrophica. Furthermore, 10 D1 and 10 D2 sequences from photosynthetic eukaryotes from taxa with a well characterized fossil record were added to allocate calibration points. The phylogeny of Type II reaction centers was crosscalibrated on D1 and D2 as listed in Table 4 and calibration points were assigned as presented in Fig.   4, red dots.
Dates for the Arabidopsis/Populus divergence, the divergence of the angiosperms (Amborella), gymnosperms (Cycas), and land plants (Marchantia) were implemented as suggested and discussed by Clarke et al. (2011), representing points 1 to 4 in Fig. 4. Three ages from eukaryotic algae were used too: an age of 190 Ma was assigned to the divergence of Phaeodactylum trichornutum and Thalassiosira pseudonana, based on fossil Jurassic diatoms from the Lias Group, reviewed by Sims et al. (2006). A minimum age of 600 Ma based on a Late Neoproterozoic Chinese multicellular red alga fossil (Xiao et al., 2004) was assigned to the split between the diatom sequences and the sequences from Porphyra purpurea, as a conservative estimate for the divergence of complex red algae, which predates this time (E. C. Yang et al., 2016). The earliest widely accepted 1.2 Ga redalgae fossil, Bangiomorpha (Butterfield, 2000;Knoll et al., 2013), was assigned as a minimum age to the divergence of the sequences from Cyanidium caldarium, a unicellular early-branching red algae, this was used as a conservative estimate for the most recent common ancestor of red algae. It is however quite likely that the last common ancestor of red algae significantly predates this age as suggested by recently described 1.6 Ga red algae fossils (Bengtson et al., 2017) and recent molecular clock analysis (Sánchez-Baracaldo et al., 2017;E. C. Yang et al., 2016).
The earliest well-assigned filamentous Cyanobacteria fossils of comparable size to those of the earlybranching Pseudanabaena have been reported at 1.9 Ga (Golubic & Lee, 1999;Sánchez-Baracaldo et al., 2017;Schirrmeister et al., 2016;Sergeev et al., 2002), and this was assigned as a minimum age to the sequences from Pseudanabaena biceps.
The oldest calibration point, point 11, was selected to be the branching point of the D2 and the G4 D1 from Gloeobacter violaceous. This was set to be around the age for the GOE and was assigned as a minimum age of 2.45 Ga (Calibration 1) (Bekker et al., 2004). For comparison, a calibration of 2.7 Ga was also used (Calibration 2) to test the effect on the estimated divergence times of an older age for crown group Cyanobacteria. Geological evidence suggests that oxygen 'whiffs' or 'oases' could significantly predating the GOE (Lyons et al., 2014) so this scenario is not entirely implausible.
The calibration points on D1 were allocated on Group 4 because this type of D1 is the only one retained in all Cyanobacteria with PSII, it is the only type of D1 inherited by photosynthetic eukaryotes, and it is the main D1 used to catalyze water oxidation under most growth conditions. The analysis by Cardona et al. (2015) indicated that the gene duplication events that led to the atypical forms and the microaerobic forms of D1 are likely to predate the last common ancestor of Cyanobacteria.
It is well accepted that a form of anoxygenic photosynthesis had already evolved by 3.5 Ga. This is demonstrated by both sedimentological and isotopic evidence for photoautotrophic microbial communities recorded in Paleoarchean rocks (Butterfield, 2015;Nisbet & Fowler, 2014;Tice & Lowe, 2004). In addition, sedimentary rocks and banded iron formations from Isua, Greenland, hint at the presence of photosynthetic bacteria in the marine photic zone as early as 3.7-3.8 Ga (Czaja et al., 2013;Grassineau et al., 2006;Knoll, 2015;Rosing, 1999;Rosing & Frei, 2004;Schidlowski, 1988). Therefore, we used a root prior of 3.5 Ga as a conservative estimate for photoautrophy based on Mgtetrapyrrole containing photochemical reaction centers. Nevertheless, because it is not yet known exactly when photosynthesis originated for the first time we also tested the effect of varying the root prior from 3.2 to 4.1 Ga on the estimated divergence time under restrictive and flexible scenarios.
A Bayesian Markov chain Monte Carlo approach was used to calculate the estimated divergence times. We used Phylobayes 3.3f , which allows for the application of a relaxed log-normal autocorrelated molecular clock under the CAT + GTR + Γ model (Lartillot et al., 2009;Lartillot & Philippe, 2004) necessary for the implementation of flexible boundaries on the calibration points (Z. Yang & Rannala, 2006). To understand the effect of different evolutionary models on the age estimates we compared the CAT + GTR + Γ model with 1) a LG + Γ model that sets less flexible boundaries on the calibration points, 2) a CAT + Γ model assuming a uniform (Poisson) distribution of amino acid equilibrium frequencies, or 3) an uncorrelated gamma model where the rates of substitution can vary independently. The flexible bounds on the CAT + GTR + Γ model were set to allow for 2.5% tail probability falling outside each calibration boundary or 5% in the case of a single minimum boundary. All molecular clocks were computed using four discrete categories for the gamma distribution and four chains were run in parallel until convergence.
In this work we define the period of time between the duplication event that led to the divergence of D1 and D2 and the appearance of the ancestral standard D1 as Δ T. This value is calculated as the subtraction of the mean age of the latter node (Fig. 4, green dot) from the former's mean node age (Fig. 4, D0, orange dot). The instant rates of evolution, which are necessary for the computation of divergence time from the phylogeny, were retrieved from the output files of Phylobayes. These rates are calculated by the software as described by the developers elsewhere (Kishino et al., 2001;Lepage et al., 2007) and are expressed as amino acid changes per site per unit of time. The rate at node D0 was termed ν max and a baseline rate of evolution during the Proterozoic was obtained as the average of all node rates in Group 4 D1 and D2 and denoted ν min . All sequence datasets and estimated divergence times for each node of each tree used in this analysis are provided in the supplementary data files.   Type II reaction center proteins. (b) An schematic representation of the phylogeny shown in (a). All reaction centers have common ancestry and descended from a homodimeric reaction center, marked A. From A, two new reaction centers emerged, one ancestral to all Type II reaction centers (II) and a second ancestral to all Type I reaction centers. This is the earliest diversification event of reaction center proteins that can be inferred from sequence and structural data and it is marked 1. The ancestral Type II reaction center protein (II) gave rise to two new proteins, one ancestral to D1 and D2, named D0 and a second ancestral to L and M named K. The ancestral L and M subunits further diversify into Chloroflexi-type and Proteobacteria-type L and M subunits (5).
Step 6 indicates that Type I reaction center proteins also diversified in parallel to Type II reaction center proteins. (c) Evolution of cyanobacterial D1 and D2, modified from Cardona et al. (2015). The conserved structural and functional traits between D1 and D2 indicate that the ancestral homodimeric photosystem was able to oxidize water, this is represented by a fuchsia spot in each D0 monomer. The lack of fuchsia dot in the G0 sequence indicates the absence of some of the ligands to the Mn 4 CaO 5 cluster. This sequence has acquired changes that are predicated to impair function and it is not known whether it is assembled into PSII. The pentagon in G1 indicates the possible binding of chlorophyll a as hypothesized by M.
Y. Ho et al. (2016). The green triangle in the G2 sequences indicates the possibility of binding of a modified cluster as these sequences still retain some of the ligands to the Mn 4 CaO 5 cluster. cluster in D1 and the corresponding site in D2. The residues colored in gray are provided by D1 and D2 protein and those colored in orange by CP43 and CP47. Notice that while a redox Y-H pair is conserved in D1 and D2, D2/CP47 evolved to stop the binding of metals by replacing ligands for phenylalanine residues creating a hydrophobic patch. Such a site is not found in any other Type II or Type I reaction center protein and it is unique to D2 . (c) Overlap of D1 and D2 highlighting some of the conserved traits as described in the text  supplementary Table S2. The gray dots around 3.5-3.8 Ga marks an approximated timing for the earliest events in the history of photosynthesis: the divergence of D1 and D2 (~29% sequence identity), the divergence of anoxygenic (L/M) and oxygenic (D1/D2) reaction center proteins (~17%), and the divergence of Type I and Type II reaction center proteins (≤5%). The curved blue line highlights that any scenario for the diversification of reaction centers after the origin of life requires faster rates of evolution at the earliest stages in the evolution of photosynthesis and D2 nodes with changed root priors, the red curve farthest to the right was calculated using a root prior of 4.2 ± 0.05 Ga, while the green curve farthest to the left was calculated using a root prior of 0.8 ± 0.05 Ga. (c) Change in the rate of evolution as a function of Δ T, the dashed lines represents a fit to a power-law function. The fit components are shown in supplementary Table S5 Figure 6. Scenarios for the evolution of Type II reaction centers. The rows of colored dots represent estimated divergence times calculated using the CAT + GTR + Γ and LG + Γ models and root priors of 3.8 or 3.5 ± 0.05 Ga as indicated in the illustration. Primordial forms of water oxidation emerged before the gene duplication event that led to D1 and D2 in an ancestral homodimeric photosystem (orange dot). The more ancient the ancestor of crown group Cyanobacteria is (green arrow), the closer the origin of water oxidation has to be, relative to the origin of the first photochemical reaction centers in the early Archaean (orange arrow). The yellow vertical bar marks the GOE