Possible link between Earth’s rotation rate and oxygenation

Modelling microbenthic O₂ export

We explored how microbial processes and export fluxes of their metabolic substrates and products from ancient benthic photosynthetic ecosystems were influenced by daylength, environmental conditions and various regulatory mechanisms of photosynthetic production and respiration using an in silico microbenthic model. Model scenarios were constructed and simulated using MicroBenthos software¹². MicroBenthos model definitions and parameters for the described scenarios are provided with this article. The software and usage instructions are available at https://microbenthos.readthedocs.io.

The modelling framework is an adaptation of de Wit et al.⁶¹. Briefly, benthic systems are constructed as a diffusive–reactive system in a 1D computational domain, with discrete cells used to represent the spatial distribution of the state and parameter variables. While the study by de Wit et al.⁶¹ focused on biomass growth running over long simulation times, our interest was to study the dynamics of process rates and solute fluxes over diel timescales. Therefore, we set a fixed biomass for the microbial groups, added a water subdomain on top of the sediment as a diffusive boundary layer and ran simulations until a diel steady state was reached (5 days). Our model domain used 5 µm cells, with an 8 mm sedimentary subdomain and 1 mm diffusive boundary layer of water on top. O₂ and sulfide concentrations were the state variables that we solved for. Photosynthetically active radiation (PAR) was expressed as a percent of the maximum intensity at the diel zenith, and followed a cosinusoidal pattern similar to that of diel insolation dynamics.

R_aero and SOX were formulated to occur throughout the sediment. Microbial groups (cyanobacteria and SRB) were represented as biomass distributions in the sediment subdomain, and biomass-dependent metabolism kinetics were expressed as multiplications of the response functions of salient environmental and state variables. Coupled partial differential equations of the state variables (O₂ and H₂S) were composed from the reaction terms that accounted for sediment porosity and were solved with finite-volume numerical approximations⁶².

Our in silico mat allowed us to explore how diffusive mass transfer shapes the interplay between illumination dynamics, gross production and consumption rates, and diel O₂ export. The effect of daylength was studied by varying the period of the illumination from 12 h to 24 h, the range of estimated daylengths over Earth’s history after the earliest estimates for the origin of OP⁶³. We report the calculated average diel net export and process rates in units of mmol m⁻² h⁻¹ because the hour is the largest temporal unit unaffected by changes in the Earth’s rotation and thus allows for comparison across daylengths.

First, we explored the simplest case of O₂ production, which is with light availability. Two microbial processes were considered: OP performed by cyanobacteria and R_aero. The parameters for the biotic reactions were re-expressed as a biomass-specific maximal yield (Q_max). A fundamental assumption is that the photosynthesis rate is strictly correlated to the instantaneous photon flux:

where sat is a Michaelis–Menten function with K_PAR = 15% and the cyanobacterial biomass with a log-normal distribution with a peak value of 12 mg cm⁻³ at 0.5 mm depth (Supplementary Video 1). The only source of O₂ is OP, and the sinks are aerobic (sedimentary) respiration (R_aero). For the production and consumption rates of C_org, we assumed a stoichiometry of:

with respect to O₂ cycling rates, where CH₂O refers to one C_org equivalent. Assuming that C_org is predominantly particulate, with negligible diffusional transport, diel C_org burial was thus calculated as:

where ∫OP and ∫R_aero are the diel depth-integrated rates of O₂ production and consumption and are equivalent to C_org production and consumption according to equation (2). Thus, diel burial can also be represented through the export flux of O₂ at the top and bottom interfaces of the sedimentary domain:

which allowed us to assess the dynamic steady state of the diel model when the average diel depth-integrated rates equalled the export fluxes.

To calibrate the O₂ productivity for unitless PAR intensities, we determined the Q_max that caused a maximum O₂ export that corresponded to the median maximal flux from illuminated benthic photosynthetic systems¹³. A Q_max of 4.0022 mmol g⁻¹ h⁻¹ produced the target export flux of 5.76 mmol m⁻² h⁻¹ under a sedimentary respiration load of 0.1 mM h⁻¹. Note that by calibrating the productivity to the maximum diel illumination, the model represents a ‘mean solar day’ of a given Earth year⁵⁹. This allowed us to disentangle the effect of daylength from geological-scale changes in the insolation intensity, such as in the ‘faint young Sun’ paradigm reviewed thoroughly by Feulner²³, or changes in the solar spectrum related to atmospheric composition⁶⁴.

Next, we explored the effect of R_anaero on the daylength dependency of the process rates and export fluxes. We used the example of sulfate reduction performed by SRB with a log-normal biomass distribution with a peak value of 2 mg cm⁻³ (Supplementary Video 1). The R_anaero rate was either defined as a constant rate process for the scenario ‘OP SRB constant’ as:

or as an O₂– and H₂S-sensitive process as:

where inhibition is a function of the local H₂S and O₂ concentration (x) of the form:

when x < K_max and 0 when x ≥ K_max. Inhibition factors chosen for both scenarios with O₂– and H₂S-sensitive SRB (‘OP SRB’ and ‘OP SRB inhibited’) were 3 mM for (K_{{rm{max}},{rm{H}}_2{rm{S}}}), 0.5 mM for (K_{{rm{half}},{rm{H}}_2{rm{S}}}) and 1 µM for (K_{{rm{half}},{rm{O}}_2}). (K_{{rm{max}},{rm{O}}_2}) was 0.8 mM for the scenario with a moderate inhibition, ‘OP SRB’, and 0.3 mM for the scenario with a strong inhibition, ‘OP SRB inhibited’.

SOX was formulated as:

where k = 351 l mol⁻¹ h⁻¹ (ref. ⁶⁵). For R_anaero we assumed the stoichiometry:

and therefore calculated diel burial as:

where ∫R_anaero is the depth-integrated rate of sulfide production by SRB. The export flux of O₂ is shaped by mat-intrinsic rates of R_aero, OP and SOX according to:

where ∫SOX is the depth-integrated rate of O₂ consumption by SOX. Assuming the complete oxidation of H₂S to sulfate, H₂S export can be formulated as:

According to equations (10)–(12), diel burial can thus be represented as:

This illustrates that the control of diel burial is related to the export of O₂ and the reduced product of R_anaero (such as H₂S), as the former is an equivalent source and the latter an equivalent sink of C_org within the mat. This means that an increase in O₂ export would not result in an increase of burial if H₂S export increased proportionally in terms of C_org equivalents.

To calibrate the productivity for SRB scenarios, we determined values of Q_max for OP and SRB, which yielded the target maximum export flux of 5.76 mmol m⁻² h⁻¹ at 24 h under 250 µM O₂ boundary conditions and negligible burial fluxes at 18 h under anoxic boundary conditions, that is, in pre-GOE conditions (see Supplementary Data 1 for model parameters).

We tested the sensitivity of diel burial to O₂ concentration in the water column (0–250 µM) for all three SRB scenarios. For the least O₂-sensitive scenario (‘OP SRB’), we also tested sensitivity of burial to gross productivity by varying the maximal photosynthetic yield Q_max over the range 0.5–10 mmol g⁻¹ h⁻¹. Note that the variation in Q_max can be considered equivalent to variations in other factors that influence gross production, such as nutrient supply and irradiance levels.

We then explored the effect of AP and reductant supply to the mat. Reductants that served as electron donors for AP were available in Precambrian phototrophic habitats, and supported diverse forms of photosynthesis even after the GOE and the evolution of OP^66,67,68. Although the extent and ecological niches of AP and OP over Earth’s history remain unclear, AP and OP probably co-existed, with spatially and temporally variable partitioning of the total GPP between them⁶⁹. Analogous to modern systems, the partitioning probably depended on the limiting factors of both metabolisms, such as light and nutrients, with AP additionally limited by electron donor supply, and on the onset of novel evolutionary avenues or geochemical transitions that facilitated shifts in the outcome of competition.

We implemented this concept by adapting the ‘OP SRB’ scenario to include metabolic flexibility in the modelled photosynthetic community, analogous to cyanobacteria that can partition harvested light energy towards OP and sulfide-driven AP (Extended Data Fig. 4). Transitions between photosynthetic modes are based on the local sulfide and light availability with the rate of OP regulated by the rate of AP²⁰. In this ‘OP SRB + AP’ scenario, the modelled cyanobacteria produced O₂ according to:

where normsat is a normalized Michaelis–Menten function with a value of 1 for H₂S ≥ H₂S_thr, where H₂S_thr is the threshold sulfide level. Conversely,

The metabolic behavior is such that OP occurs only below H₂S_thr, which is light dependent. Below this threshold level, the harvested light energy is partitioned towards AP and OP, with a higher affinity for AP. Above the threshold, only AP occurs and OP is suppressed. The resulting metabolic response is that the sum of OP and AP follows the form of equation (1), whereas R_anaero, R_aero, SOX and other processes work as in previous scenarios. The interplay of illumination, mat processes and the resultant depth-resolved dynamics of O₂ and H₂S under this scenario can be seen in the Supplementary Video 1. Assuming sulfate as the product of AP according to:

diel burial in this scenario was calculated as

where ∫AP is the depth-integrated rate of CO₂ fixation by AP. As in the other scenarios, diel burial can alternatively be calculated according to equation (13). Values of Q_max for OP and R_anaero were adjusted to produce the same maximum O₂ export flux at 24 h (O₂ boundary at 250 µM) and diel C_org flux at 21 h (25 µM O₂ boundary) as in ‘OP SRB’. This allowed us to establish directly comparable upscaling calculations for global (p{{rm{O}}_2}) at the same areal coverage and diel burial by ‘OP SRB’ and ‘OP SRB + AP’ mats (see below).

Besides the effect of daylength, we also studied the effect of the availability of O₂ and H₂S as electron donors from the water column on the scenario with AP. Water column O₂ levels from 0 to 250 µM and H₂S levels from 0 to 20 µM were tested.

Photosynthetic inhibitory mechanisms were also tested by modelling cyanobacteria that exhibit a 30 min recovery time of OP after the local sulfide levels fall below a certain level²⁷.

Microsensor measurements

We sampled cores from cyanobacterial mats that formed in the MIS (Michigan, United States) in October 2015, May 2016 and June 2016. Cores and bottom water sampled from above the mats were transported to the lab in Ann Arbor. Cores were either used directly or kept on the seasonal day–night light cycle at 8 °C until the measurements were made. During the measurements in the cores, a circular flow of the water column above the mat surface was adjusted using peristaltic pumps connected to an external temperature-controlled reservoir of bottom water. Water column O₂ and sulfide concentrations were adjusted in this separate reservoir of bottom water using N₂, air and neutralized Na₂S solution. After adjustment, the reservoir was covered with paraffin oil to prevent gas exchange and maintain constant conditions in the water column. The total sulfide concentration in the water column was regularly checked by subsampling and colorimetric determination according to Cline⁷⁰. Both the water column reservoir and the core were kept at 10 °C during the measurements.

Illumination was provided by a broadband halogen lamp (Schott KL-1500). In situ spectral light quality was simulated using optical filter foils (Roscolux 375, Rosco Laboratories). The incident irradiance at the surface of the mat was determined with a submerged cosine-corrected quantum sensor connected to a LI-250A light meter (both LI-COR Biosciences GmbH). O₂, pH and H₂S microsensors with a tip diameter of 10–50 µm and a response time of <1 s were built, calibrated and employed as previously described^71,72,73. The microsensor tips were always separated by <500 µm during simultaneous O₂, pH and H₂S measurements.

The volumetric rates of gross OP were estimated based on the dynamics of O₂ concentration after a light–dark shift, as described previously⁷⁴. Analogously, the volumetric rates of gross AP were calculated from the increase of H₂S concentration and pH directly after a light–dark transition (that is, light-driven S_tot consumption rates), as previously described^26,27,75. Net rates and fluxes were calculated from profiles based on Fick’s laws using diffusion coefficients for O₂ and sulfide corrected for temperature and salinity.

From diel export to global oxygenation

We considered that changes in diel O₂ and H₂S export rates, the difference of which is equivalent to diel benthic C_org burial, interact with global redox controls beyond the mat domain (Extended Data Fig. 1). Namely, we considered atmospheric reduction by volcanism- and metamorphism-derived gases (atmR) and weathering⁶⁰. Although we assumed a constant flux of reductant (vR) available for atmospheric reduction, weathering is dependent on (p{{rm{O}}_2}) and thus strongly tempers the accumulation of O₂ in the atmosphere. As a recent computational Earth rotation model¹⁰ predicts a stabilized daylength in the mid-Proterozoic, we started our global calculations by assuming a (p{{rm{O}}_2}) of 0.01–0.1 at the predicted stable 21 h. This (p{{rm{O}}_2}) is the result of the interaction of global sources and sinks (illustrated in Extended Data Fig. 1) according to a simplified version of the formulations (for example, COPSE^60,76) of the global O₂ budget. Considering an evolution of steady states of (p{{rm{O}}_2}) throughout Earth’s history, we described the (p{{rm{O}}_2}) evolution at a given age as:

which provides the steady state (p{{rm{O}}_2}) as:

where pB is the global burial flux associated with marine pelagic production, mB is total mat burial flux, vR is the flux of volcanic reductant, tB is the burial flux by terrestrial mats, and uB is an aggregate flux term that captures uplift forcing, the global C_org reservoir and a weathering constant⁶⁰.

Changes in any of the global fluxes would therefore result in new steady state (p{{rm{O}}_2}). For the total marine burial (pelagic pB and coastal benthic mB – tB) we assumed modern values throughout the Precambrian⁶⁰. Beyond the marine realm, we considered the possible existence of terrestrial mats⁴³, which would have further boosted global GPP and C_org burial, but would have also been more susceptible to weathering (tB in Extended Data Fig. 1). For these terrestrial mats, we assumed the same regulation mechanisms of export fluxes as those for coastal benthic mats and considered that AP might have been driven by reductant supplied from R_anaero, but we did not explore the effects of external reductant availability for terrestrial mats. Using our diel mat model output for tB at 21 h daylength and for 5% continental coverage, we then calculated vR for the two (p{{rm{O}}_2}) levels assuming the values for uB from Daines et al.⁶⁰.

The Earth rotation rate model¹⁰ predicts monotonic deceleration before 2.2 Ga followed by a stable daylength of ~21 h due to the atmospheric thermal tide resonance until ~650 Ma, with a subsequent increase towards the modern 24 h daylength¹⁰. Starting at 21 h at 2 Gyr, we estimated the effect of the daylength-driven change in the C_org burial from benthic and terrestrial mats, the output of our diel mat model, on (p{{rm{O}}_2}) levels both backwards and forwards in time. The modulation of (p{{rm{O}}_2}) naturally depends on the steepness of the dependence of mat C_org burial on daylength (Fig. 2). However, we limited our analysis to the effect of mats with metabolic regulation from the moderately daylength-sensitive scenarios ‘OP SRB’ and ‘OP SRB + AP’. Scenarios with a steeper dependency of burial on daylength yielded oscillations of (p{{rm{O}}_2}) beyond PAL even due to the minimal variations in daylength within the resonant-locked phase in the mid-Proterozoic¹⁰ and therefore had to be excluded from analysis of (p{{rm{O}}_2}) evolution. We additionally included negative feedback effects of increasing (p{{rm{O}}_2}) on the mat C_org burial by dynamically adjusting the O₂ boundary conditions in each time step according to the (p{{rm{O}}_2}) calculated from the previous step. To then calculate the quasisteady-state (p{{rm{O}}_2}) level based on changes in the benthic and terrestrial burial, the actual coverage and thus partitioning between the pelagic and benthic burial have to be considered. Reliable estimates for this partitioning are lacking. The presence of pelagic cyanobacteria has persisted since 1.1 Ga (ref. ³⁸). Yet, cyanobacteria in the palaeontological record are primarily benthic, with sparse evidence for pelagic forms^39,40,41. Hypotheses for the limited pelagic productivity during the Proterozoic range from an inaccessible, toxic or damaging photic zone to late evolution of a planktonic lifestyle^21,66,77. Alternatively, burial of the C_org produced by pelagic cyanobacteria might have been hindered by the small cell sizes, low sinking rates and resultant efficient remineralization within the water column^42,78. The preservation of C_org by benthic microbial mats is as uncertain, especially when eroded—on land or in the oceans. In our calculations we considered the terrestrial erosional weathering explicitly (tB in Extended Data Fig. 1), which introduces a strong negative feedback effect on (p{{rm{O}}_2}). This implicitly describes a substantial loss in translation of the diel C_org to long-term burial that we, however, did not explicitly account for. For marine benthic mats we argue that a more direct link between the diel and long-term burial of C_org is plausible given that mats can reach substantial thicknesses when undisturbed (for example, in Solar Lake the thickness is >1 m) (ref. ⁷⁹) and that mats have a similar remineralization fate as pelagic C_org export when eroded. To address these uncertainties, we considered benthic burial as between 20 and 50% of the total marine burial at 21 h daylength. Note that this partitioning must shift during our simulations over Earth history because benthic (and terrestrial) burial is daylength dependent, whereas we took an ‘all is constant’ approach for pelagic burial and the fraction of weathered diel C_org, as we expected a negligible effect of molecular diffusion and thus daylength on these fluxes.

Source: Ecology - nature.com