The Great Oxygenation Event as a consequence of ecological dynamics modulated by planetary change

The Great Oxygenation Event (GOE), ca. 2.4 billion years ago, transformed life and environments on Earth. Its causes, however, are debated. We mathematically analyze the GOE in terms of ecological dynamics coupled with a changing Earth. Anoxygenic photosynthetic bacteria initially dominate over cyanobacteria, but their success depends on the availability of suitable electron donors that are vulnerable to oxidation. The GOE is triggered when the difference between the influxes of relevant reductants and phosphate falls below a critical value that is an increasing function of the reproductive rate of cyanobacteria. The transition can be either gradual and reversible or sudden and irreversible, depending on sources and sinks of oxygen. Increasing sources and decreasing sinks of oxygen can also trigger the GOE, but this possibility depends strongly on migration of cyanobacteria from privileged sites. Our model links ecological dynamics to planetary change, with geophysical evolution determining the relevant time scales.

T he Great Oxygenation Event (GOE),~2.4 billion years ago, records a major turning point in the history of our planet. While pO 2 may have fluctuated during the GOE 1,2 , its results are clear: multiple lines of geological and geochemical evidence document the initial rise of O 2 to permanent prominence in the atmosphere and surface ocean 3 . These include (1) a sharp drop in iron formation deposition 4 , (2) the appearance of red beds in continental sedimentary successions 5 and calcium sulfates in marine environments 6 , (3) the retention of iron in ancient weathering horizons 7 , (4) the loss of detrital uraninite and other redox-sensitive minerals from fluvial and deltaic sandstones 8 , and (5) the loss of a mass-independent sulfur-isotopic signature best explained in terms of photochemical reactions in an essentially oxygen-free atmosphere 9,10 . Detrital pyrite and uraninite suggest that, prior to the GOE, pO 2 was <1.6 × 10 −4 of present atmospheric levels (PAL) 11 , whereas modeling of mass-independent S-isotopic fractionation limits Archean pO 2 tõ 10 −5 of PAL or lower 9,12,13 . Figure 1 summarizes the main geological and geochemical lines of evidence for the GOE.
While there is consensus on when the biosphere began to accumulate oxygen, debate continues about the physical and/or biological drivers of this transition. As emphasized by Kasting 14 , the GOE required an increase in the rate of oxygen production and/or a decrease in rates of oxygen consumption. Cyanobacterial photosynthesis is generally accepted as the key source of oxygen, but many models to explain the GOE tacitly assume that cyanobacteria were abundant prior to the GOE and so rely on physical events as the proximal drivers of environmental state change [14][15][16] . Proposed physical drivers include hydrogen escape to space associated with the photodissociation of H 2 O and CH 4 in the upper atmosphere 17,18 , a temporal shift to more oxidized volcanic gases as subaerial volcanism increased along with craton expansion 19,20 , or a decrease in hydrothermal iron fluxes into the ocean as vent temperatures declined 19,21,22 . Limited phosphorus availability in Archean oceans has also been suggested as limiting rates of photosynthesis and carbon burial [23][24][25] .
Goldblatt et at. 26 considered that oxygenic photosynthesis was well established prior to the GOE but that atmospheric methane oxidation suppressed oxygen levels. They modeled a low-level steady state of oxygen and a high-level steady state-the latter being characterized by an ozone layer that shields the troposphere from ultraviolet radiation, limiting the rate of methane oxidation. These steady states of oxygen can overlap, resulting in bistability and hysteresis. Also related to methane oxidation, Konhauser et al. 27 proposed a coupled biological-geological driver for the GOE, concluding that a reduced flux of nickel to the oceans in the late Archean limited methanogen activity, thereby capping the supply of biogenic methane.
Other investigations of the carbon isotope record have argued that the δ 13 C value of marine carbonate increased from the Archean into the Proterozoic 28,29 , raising the possibility that biological innovations acted as triggers for the GOE. In contrast to geophysical models, Ward et al. 30 modeled the GOE as a direct response to the origin of oxygenic photosynthesis. This implicitly assumes that cyanobacteria hold an inherent selective advantage over anoxygenic photoautotrophs, but the dominance of  anoxygenic photosynthetic bacteria in modern environments where light is present but oxygen is not (e.g., ref. 31 ) casts substantial doubt on this assumption. The competitive advantage of anoxygenic photosynthetic bacteria over cyanobacteria may relate to the lower metabolic cost of deriving electrons from donors other than water, or it may reflect the direct inhibition of oxygenic photosynthesis by sulfide 32 or ferrous iron 33 . In the presence of sulfide, some cyanobacteria can shut down Photosystem II and use H 2 S as an electron donor for anoxygenic photosynthesis 34,35 . The inhibitory effect of Fe 2+ is less clear, as Ward et al. 36 show that cyanobacteria live in hot spring waters that contain both ferrous iron and modest amounts of oxygen. Where Fe 2+ was present and O 2 was absent, however, cyanobacteria were not ecologically important 36 . Moreover, a growing array of geochemical data is interpreted as evidence for "whiffs of oxygen"-spatially and temporally limited oxygen oases in the Archean biosphere-requiring oxygenic photosynthesis hundreds of millions of years before the GOE [37][38][39] . A natural scenario, then, is that cyanobacteria were present but relatively scarce well before the GOE, rising to global ecological prominence coincident with the environmental transformation. Jones et al. 40 argued that as the ratio of alternative electron donors (they focused on Fe 2+ , but it is the sum of non-H 2 O electron donors that is key) to phosphorus declined through time, environmental opportunities for oxygenic photoautotrophs increased. Knoll and Nowak 41 developed a simple mathematical model to describe the competition between anoxygenic photosynthetic bacteria and cyanobacteria on a changing Earth. Their model keeps track of the abundances of anoxygenic photosynthetic bacteria, cyanobacteria, ferrous iron, and oxygen. Simulations revealed that the GOE can be triggered as the planetary removal rate of oxygen from the atmosphere declines, although that primitive model was not analytically solved. Ozaki et al. 42 used both modeling and experimental evidence, arguing once again that oxygenic photoautotrophs evolved long before the GOE, but that geochemical conditions throughout the Archean favored primary production by anoxygenic photosynthetic bacteria. They further argued that anoxygenic photoferrotrophs, being well adapted to low-light conditions, inhabited the bottom of the photic zone, thereby diminishing the supply of upwelling nutrients to cyanobacteria (see also ref. 43 ). This control on the proliferation of oxygenic bacteria would have remained in place as long as there was a sufficiently large supply of reducing agents from the deep ocean.
Here, we advance the model of Knoll and Nowak 41 by additionally keeping track of the abundance of phosphate, because the amount of ferrous iron or other alternative electron donors relative to phosphate influences the competition between anoxygenic photosynthetic bacteria and cyanobacteria. We also properly account for the loss of ferrous iron due to the proliferation of anoxygenic photoautotrophs, and we correctly incorporate loss of dioxygen due to iron oxidation. We can now investigate the effects of a declining influx of iron and an increasing influx of phosphate as consequences of planetary change, while concurrently probing the influences of biological innovations and of changes in sources and sinks of oxygen. Our model thus focuses consideration of the GOE on interactions between the physical and biological Earth.

Results
Based on the present-day distribution of photosynthetic bacteria 31 , we assume a competitive advantage for anoxygenic photosynthetic bacteria in early environments where electron donors such as Fe 2+ , H 2 S, or H 2 were present. We also assume the contemporaneous existence of environments where cyanobacterial populations could thrive, providing a seedbed for migration. Non-marine waters provide an example of the latter, supported by the branching of non-marine taxa from basal nodes in cyanobacterial phylogenies 44,45 and also by the presence of stromatolites in Archean lacustrine successions 46 , despite the likelihood that many Archean lakes and rivers had low levels of potential electron donors such as Fe 2+ and H 2 S 47 .
Following Jones et al. 40 and Ozaki et al. 42 , we use Fe (iron) and P (phosphorus) to represent the environment, which is similar to the H 2 and P employed in other studies 48,49 . The logic of this choice is that in Archean oceans, Fe 2+ is thought to have been the principal electron donor for anoxygenic photosynthesis 50,51 , whereas P governed total rates of photosynthesis. (Kasting 14 argued that H 2 was key to photosynthesis on the early Earth, a view supported by low iron concentrations in some early Archean stromatolites 52 .). In any event, under the conditions of low P availability thought to have characterized early oceans 25,40,49,[53][54][55] , anoxygenic photosynthesis would have depleted limiting nutrients before alternative electron donors were exhausted. In consequence, rates of photosynthetic oxygen production would be low. As iron availability declined and/or P availability increased, the biosphere would inevitably reach a point where P would remain after Fe 2+ had been depleted, expanding the range of environments where cyanobacteria are favored by natural selection 42 .
Our model keeps track of the abundances of anoxygenic photosynthetic bacteria (APB), x 1 , cyanobacteria, x 2 , and three crucial chemicals: iron(II) (Fe 2+ ), y 1 , phosphate (PO 4 3− ), y 2 , and dioxygen (O 2 ), z. Both types of bacteria require phosphate for reproduction. APB needs iron(II) (or some other suitable reductant) as an electron donor in photosynthesis. The following five equations describe the reproduction and death of APB and of cyanobacteria as well as the dynamics of iron(II), phosphate, and dioxygen: Here, we have omitted to write symbols for those rate constants that, for understanding the GOE, can be set to one without loss of generality (Supplementary Note 1). Each remaining rate constant is a free parameter. Equations (1) thus satisfy redox balance by construction. We are left with a system that has five main parameters: c specifies the rate of reproduction of cyanobacteria; f 1 and f 2 denote the rates of supply of iron(II) and phosphate, respectively; a denotes biogenic production of oxygen; b denotes geochemical consumption of oxygen. Note that iron(II) and phosphate are also removed by geochemical processes at a rate proportional to their abundance. In addition, iron(II) is used up during anoxygenic photosynthesis, and iron(II) reacts with oxygen and is thereby removed from the system. Phosphate is used up during the growth of APB and cyanobacteria. (We investigate extensions of the model that incorporate bounded bacterial growth rates and organic carbon in Supplementary Note 2 and Supplementary Note 3, respectively.) We posit iron(II) as the primary electron donor for anoxygenic photosynthesis, and for simplicity of presentation, we refer to y 1 and f 1 in this context. However, as noted above, y 1 and f 1 can similarly represent the abundances and influxes of other alternative electron donors, especially dihydrogen (H 2 ) 56,57 and hydrogen sulfide (H 2 S) 58 . Our model, its analytical solution, and the conclusions that follow hold equally well by considering any of these electron donors or all together.
We also include small migration rates, u 1 and u 2 , which allow for the possibility that APB and cyanobacteria persist in privileged sites from which they can migrate into the main arena of competition. On the Archean Earth, these parameters could have been affected by the flow of water and by surface winds. For the mathematical analysis presented in the main text, we assume that these rates are negligibly small.
The GOE represents the transition from a world dominated by APB (Equilibrium E 1 ) to one that is dominated by cyanobacteria (Equilibrium E 2 ) (Figs. S1, S2). On a slowly changing planet, the abundances of APB and cyanobacteria and of the three chemicals are approximately in steady state. Therefore, we consider the fixed points of Eqs. (1).
Pure equilibria. In the absence of APB and cyanobacteria, the abiotic equilibrium abundances of iron(II) and of phosphate are given by f 1 and f 2 , respectively, and there is no oxygen in the system. If f 1 f 2 > 1, then APB can emerge. Subsequently, the system settles to Equilibrium E 1 , where only APB are present and there is still no oxygen. E 1 is stable against invasion of cyanobacteria if This condition can be fulfilled if the influx of iron, f 1 , is large enough, or if the influx of phosphate, f 2 , is small enough. The term on the right-hand side of the inequality is an increasing function of the reproductive rate, c, of cyanobacteria. If cf 2 > 1, then the system admits another equilibrium, E 2 , where only cyanobacteria are present and oxygen is abundant. Equilibrium E 2 is stable against invasion of APB if The left-hand side of the inequality is positive. If the right-hand side is negative (that is, if f 1 < c), then the condition certainly holds. If the right-hand side is positive, then the condition can be fulfilled if the influx of phosphate, f 2 , is large enough, or if the production of oxygen, a, is large enough. In other words, the dominance of cyanobacteria after the GOE can be guaranteed by a sufficiently large supply of phosphate or sufficiently large production of oxygen. It may or may not be possible for the proportional removal rate of oxygen, b, to become small enough for the condition to be fulfilled.
Mixed equilibrium. If Conditions (2) and (3) are either both satisfied or both not satisfied, then the system also admits an interior equilibrium,Ê. If Conditions (2) and (3) Condition (4) is understood as follows. If b is sufficiently large, then there is not enough atmospheric oxygen for rusting to render E 2 stable against invasion of APB before E 1 loses stability; the result is stable coexistence. But if b is sufficiently small, then rusting causes E 2 to become stable before E 1 becomes unstable. The critical value of b therefore depends on the input of atmospheric oxygen for Equilibrium E 2 ; it is an increasing function of the reproductive rate of cyanobacteria and of their rate of production of oxygen. If a < 1, then bistability is not possible. In this case, for Equilibrium E 2 , dioxygen is depleted by rusting before there is any significant loss of iron(II). As a result, E 2 cannot gain stability before E 1 loses stability, regardless of the values of b or c. Figure 2 shows, for different values of b, the behavior of the system as a function of f 1 and f 2 .
Transition from Equilibrium E 1 to Equilibrium E 2 . The transition between Equilibria E 1 and E 2 can be achieved by reducing the supply of iron(II), f 1 , since such a reductant is required for anoxygenic photosynthesis. When this happens, we lose the stability of E 1 and gain the stability of E 2 .
The transition is gradual if b > c(a − 1). Figure 3 shows gradual oxygenation due to decreasing f 1 . In this case, the transition occurs via the mixed equilibrium,Ê, where both types of bacteria coexist (Fig. 4). A subsequent increase in f 1 can cause APB to regain dominance (Fig. S3a).
Alternatively, if b < c(a − 1), then the transition is sudden (i.e., discontinuous). Figure 5 shows rapid oxygenation due to decreasing f 1 . In this case, E 2 is already stable before E 1 loses stability (Fig. 6). This results in bistability and hysteresis: Once the world is dominated by cyanobacteria, moderate fluctuations in the supply rate of iron would no longer change the status quo (Fig. S3b).
The effects of increasing f 2 are nearly identical to those of decreasing f 1 . Increasing the supply of phosphate results in loss of stability of E 1 and gain of stability of E 2 . This is because as f 2 rises, APB proliferate, inducing a concomitant depletion of iron reserves. The transition is gradual if b > c(a − 1) (Fig. S4) or sudden if b < c(a − 1) (Fig. S5). The critical values of f 1 and f 2 are robust to changes in u 2 (Figs. S6a, S6b).
Another possibility is that the GOE resulted from an increase in the parameter c, which denotes the reproductive rate of cyanobacteria, as affected by biological mutations. We cannot exclude the possibility that cyanobacterial performance and, therefore, primary production increased as a function of genetic innovations; however, the observation that even today oxygenic photosynthesis by cyanobacteria is limited when alternative electron donors are present places limits on such speculation. The parameter c could also be affected by geophysical or geochemical properties unrelated to oxygen consumption and independent of iron(II) or phosphate flux, such as temperature, pH, salinity, or availability of trace nutrients or other resources. The transition can be gradual (Fig. S7) or sudden (Fig. S8), depending on whether b > c(a − 1) or b < c(a − 1) when c is such that Equilibrium E 1 becomes unstable. Similar to the critical values of f 1 and f 2 , the critical value of c for triggering a GOE is robust to changes in u 2 (Fig. S6c).
Yet another possibility is that the GOE was triggered by an increase in parameter a, which measures the production rate of oxygen (Fig. S9), or by a reduction in parameter b, which denotes the proportional consumption rate of oxygen (Fig. S10). For this transition to occur, however, it is essential that u 2 is sufficiently large. Moreover, the critical values of a and b are strongly dependent on the magnitude of u 2 . If a is not large enough, then it is not possible for a reduction in b to trigger a GOE, regardless of how small b becomes (Fig. S11).
A GOE resulting from an increase in a or a decrease in b is necessarily sudden. This is because as a rises or b declines, Equilibrium E 2 , which is characterized by abundance of oxygen, may eventually gain stability, while Equilibrium E 1 remains stable. If a becomes sufficiently large or b becomes sufficiently small, then E 1 may cease to exist, and a saddle-node bifurcation results in rapid oxygenation.
Effects of migration rate. The migration rates, u 1 and u 2 , have negligible effects on the abundance of APB for Equilibrium E 1 and on the abundances of cyanobacteria and oxygen for Equilibrium E 2 . The principal effects of the migration rates are to determine the abundance of APB for Equilibrium E 2 and the abundances of cyanobacteria and oxygen for Equilibrium E 1 . As such, u 1 and u 2 control the magnitude of the decline in APB across the GOE and the magnitude of the rise in cyanobacteria and oxygen across the GOE (Figs. S12, S13).

Discussion
Our analytical investigation of ecological dynamics indicates that a switch in ecological dominance from APB to cyanobacteria Fig. 3 The GOE can be triggered by a decline in the influx of iron(II) and is gradual if b > c(a − 1). Equilibrium E 1 (APB dominate) loses stability and Equilibrium E 2 (cyanobacteria dominate) gains stability when f 1 drops below f Ã 1 and f 0 1 , respectively. We set f 2 = 80, c = 10, a = 10, b = 100, and u 1 = u 2 = 10 −3 . a We simulate Eqs. (8) from Supplementary Note 1 with α 1 = α 2 = β 1 = β 2 = 1, and we set f 1 = 100 − 40(t/10 5 ). t * denotes the time at which Equilibrium E 1 loses stability. b There is stable coexistence of both types of bacteria for f 0 would have been sufficient to spawn the GOE. Accordingly, the competitive advantage of cyanobacteria, c, the influx of suitable reductants, f 1 , and the influx of phosphate, f 2 , that appear as parameters in Condition (2) are key considerations for determining when the GOE began. As extant cyanobacteria display low fitness in sunlit but anoxic environments, this observation must condition any proposed mechanisms by which c could have substantially increased around the time of the GOE. In contrast, a decrease in f 1 and/or an increase in f 2 comprise robust mechanisms for initiating the GOE. Decreasing Fe and increasing P fluxes to the oceans are both predicted by secular cooling of the mantle (and hydrothermal systems), continental emergence, and increasing oxidant supply as the GOE began 22,24,[53][54][55] . Our analysis reveals the time at which the GOE began to be determined by the difference, f 1 − f 2 , in these influxes. It does not depend on these influxes individually. This functional dependence is preserved under the assumption that bacterial growth rates are limited (Supplementary Note 2, Fig. S14), as are the possibilities of both gradual and sudden transitions (Figs. S15, S16). These results also hold when explicitly accounting for organic carbon (Supplementary Note 3, Fig. S17). While prior investigations have focused heavily on sources and sinks of oxygen as potential drivers of the GOE, our analysis emphasizes that this possibility hinges critically on competition from cyanobacteria. If u 2 was small, then sources and sinks of oxygen, a and b, were mostly irrelevant for initiating the GOE. There is a simple intuition behind this observation. If cyanobacteria were ecologically subordinate to APB and scarce in the Archean, then atmospheric levels of O 2 would have remained low. a might have increased and b might have decreased-and fractional changes in these parameters could have been substantial-but when multiplied by low abundance of cyanobacteria, the absolute change in atmospheric O 2 levels would also be small.
Our model of ecological dynamics is robust to a broad range of influences on primary production, oxygen generation, and oxygen consumption. Our study also emphasizes that it was not strictly geophysical processes or biological innovations that ushered in the GOE, but rather the interplay between Earth and life as populations adapted to a changing planet.

Methods
We elaborate our mathematical model of ecological dynamics, its analytical solutions, the stability properties of its fixed points, and important considerations behind the GOE in Supplementary Note 1. Here, we provide an abbreviated account of our analysis and key findings.
Ecological dynamics and fixed points. Equations (1) specify the ecological dynamics of APB and cyanobacteria. They set the foundation for our understanding of the GOE. To make progress analytically, we make two simplifying assumptions. First, we assume that, at any given time, Equations (1) are approximately in steady state. This is because f 1 , f 2 , c, a, and b change very slowly relative to the typical reproductive lifetimes of the bacteria. Second, we assume that u 1 and u 2 are both small. Equations (1) become Equations (5) admit a solution for which the equilibrium abundance of cyanobacteria is zero: 2 ¼ 0 Equations (5) admit another solution for which the equilibrium abundance of APB  9 (e), and 4 (f), the stable equilibrium (green dot) moves continuously from a world that is dominated by APB to one that is dominated by cyanobacteria. Parameter values are f 2 = 10, c = 1, a = 10, b = 12, u 1 = u 2 = 1, and α 1 = α 2 = 1. The GOE is gradual.
Timing and nature of the GOE. The GOE corresponds to a transition between Equilibrium E 1 and Equilibrium E 2 . One possibility is that the GOE is gradual. Initially, E 1 is stable, while E 2 is unstable. When E 1 loses dynamical stability, a stable interior fixed point, given by Eqs. (8), (9), (10), and (11), appears near E 1 in phase space. As parameter values become more favorable to cyanobacteria, the interior fixed point moves toward E 2 in phase space, and oxygenation is progressive. When E 2 gains dynamical stability, the GOE is complete.
Another possibility is that the GOE is sudden. In this case, E 2 gains dynamical stability first. An unstable interior fixed point, given by Eqs. (8), (9), (10), and (11), appears near E 2 in phase space, and the interior fixed point moves toward E 1 in phase space. When E 1 loses dynamical stability, sudden oxygenation results, and the GOE is complete.
Numerical integration. We used the fourth-order Runge-Kutta method to numerically integrate our differential equations.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.