What Fractionates Oxygen Isotopes During Respiration? Insights from Multiple Isotopologues and Theory

The precise mass dependence of respiratory O2 consumption underpins the “oxygen triple-isotope” approach to quantifying gross primary productivity in modern and ancient environments. Yet, the physical-chemical origins of the key 18O/16O and 17O/16O covariations observed during respiration have not been tied to theory; thus the approach remains empirical. First-principles calculations on enzyme active-site models suggest that changes in the O-O bond strength upon electron transfer strongly influence respiratory isotopic fractionation. This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. 2 However, molecular diffusion may also be important. Here, we use measurements of the relative abundances of rare isotopologues 17O18O and 18O18O as additional tracers of mass dependence during dark respiration experiments of lacustrine water. We then compare the experimental results to first-principles calculations of O2 interacting with heme-oxidase analogues. We find a significantly steeper mass dependence, supported by theory, than has been previously observed. Enrichments of 17O18O and 18O18O in the O2 residue suggest that q values are strongly influenced by chemical processes, rather than being dominated by physical processes (i.e. by bond alteration rather than diffusion). In contrast, earlier data are inconsistent with theory, implying that analytical artifacts may have biased those results. Implications for quantifying primary productivity are discussed. Table of


Introduction
The isotopic composition of molecular oxygen in Earth's atmosphere is determined by chemical reactions associated with biological and photochemical O 2 cycling. 1 Molecular oxygen

Isotopic notation
We will provide a brief review of stable isotope notation before detailing the relationships between two isotopes or isotopologues during fractionation. Because the fractionation of oxygen isotopes during respiration is the primary focus, we will use notation specific to this topic as examples of these broader concepts. More general examples of these equations can be found elsewhere. 32,33 First, R describes the ratio of minor isotopes (i.e. 18 O or 17 O) to a more abundant isotope ( 16 O), and d-notation quantifies the difference of these ratios between an analyte sample and standard, reported in per mil (‰). (1) Herein, we use atmospheric O 2 as the primary standard. During a process such as respiration, the fractionation of isotopes between two phases, materials or reservoirs (i.e., respired oxygen and the remaining residue of oxygen) is described using a notation. (2) We note here that in some biogeochemical studies, e notation has also been used to describe isotope discrimination during fractionation, i.e., . 10 During a mass-dependent process, isotopic species fractionate in a co-varying manner related by the quantity θ.
This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.

(4)
It is often beneficial to linearize these equations to simplify mathematical operations, so the following expressions can be used in lieu of Equations 1-3: The linearized notation is distinguished by the inclusion of a prime (′) symbol. Expressions for q can thus be determined as (8) Finally, one may be interested in relating a substance's isotopic composition to an expected composition predicted by a theoretical or empirical mass-dependent relationship. Here, we use "l" in place of "θ" to highlight that this value is often an empirically derived relationship (i.e. the terrestrial fractionation line or closed-system respiration). This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. 7 The measure of the substance's departure from expectation, D ¢17 O (reported in parts per million, ppm), has been applied to many systems both terrestrial and planetary, and we note that other communities may refer to this parameter as " 17 O-excess" or the entire system as "triple oxygen isotopes." 15, [34][35][36][37][38] The value λ 17/18 = 0.518 is typically used for marine GOP studies. 21,39 θ values in singly substituted isotopologues can be variable Mass-dependent processes typically fractionate δ′ 17 O values approximately half as much as δ ′18 O values; therefore, θ is ~0.5. Subtle deviations from this nominal value may be used as a process diagnostic. For example, the rates of kinetic processes may scale with atomic or molecular velocity (e.g., particle mass for diffusion-limited processes). Hence, the θ values of certain kinetic processes can be approximated using the atomic or molecular mass of oxygen in Eq. 10: (10) where m is an atomic or molecular mass and subscripts refer to the three stable oxygen isotopes ( 16 O, 17 O). Solving Eq. 10 using these respective mass values yields θ 17/18, kinetic = 0.516 and θ 33/34, kinetic = 0.509 for atomic and molecular oxygen, respectively (subscripts indicate the cardinal mass of the species in the numerator and denominator of Eq. 9, respectively).
In contrast, equilibrium processes depend on the thermodynamics of bonds in the coexisting species. 40,41 The θ values associated with equilibrium processes are typically temperaturedependent, but at the high-temperature limit they approach the following value: This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. (11) .
The expression above was originally derived assuming that vibrations are harmonic, but it has nevertheless proven to be a good guidepost in a variety of systems. 32,33 Solving Equation 11 using the exact masses for atomic oxygen yields θ = 0.5305.
Specific processes differ subtly in their fractionation tendencies, however, complicating the simple view presented above. 32,33,42,43 For example, the θ value most frequently associated with dark respiration lies between 0.514 and 0.516. 12,44 This constant was empirically determined for a range of organisms and found to be nearly invariant despite much larger variability in 18 This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.

Accurate θ 17/18 values are crucial for determinations of GOP
The precise value of θ 17/18 (and λ 17/18 ) for respiration in natural systems is one of three elements that must be determined independently for triple-isotope-based estimates of GOP. The others are the isotopic composition of atmospheric O 2 dissolved in water at equilibrium and the isotopic fractionation for photosynthetic O 2 production from water. Together, these elements allow one to construct a mass balance of O 2 within the oceanic mixed layer. A brief description of the method follows.
Atmospheric O 2 has a unique ∆′ 17 O value inherited from mass-independent isotope exchange during photochemical reactions between CO 2 , O 3 and O 2 in the stratosphere. 15,45,46 This value is distinct from that produced by photosynthesis, the other main source of O 2 in the surface ocean. 3,4 At its simplest, dissolved O 2 is the combination of these two sources, and ∆′ 17 O can be used as a tracer of the mass balance between them. Respiration and gas exchange with the atmosphere are additional processes that influence ∆′ 17 O and must be accounted for. Using a λ 17/18 value of 0.518 preserves the D ¢17 O of O 2 during respiration in principle, while gas exchange tends to drive the dissolved O 2 composition toward solubility equilibrium (Figure 1a). 21 The method has been instrumental in assessing open-ocean GOP over the past two decades, offering a valuable spatiotemporal bridge between traditional short-term incubation methods and regional-to-global satellite estimates derived from ocean color. 16,18,[47][48][49] Each method has its unique uncertainties. 19, 50-52 Yet as the triple-isotope field has matured, it has become clear that many factors can exert a significant influence on D ¢17 O-based estimates of productivity, e.g., uncertainties in gas exchange rates and isotopic perturbations associated with non-steady-state conditions variable expression of respiration mechanisms and their fractionation factors; the This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
source-water isotopic composition; and mathematical approximations used in calculations. 3,12,22,26,29,[53][54][55][56][57][58][59][60] To the best of our knowledge, however, the origin of the λ 17/18 value of 0.518 has not been studied. It is consistent with the accepted range of θ 17/18 values measured for dark respiration (0.514 -0.516) and lies close to that from the simplified calculation for kinetic fractionation (see Eq. 10). 12,23,44 Yet a clear mechanism that explains this trend, rooted in physical and/or biochemical processes, has been lacking. The potential interplay of diffusion, O 2 -enzyme interactions, and O 2 reduction kinetics in the expression of isotopic fractionation suggests that multiple factors could be important for determining the respiratory θ 17/18 value. Without a mechanistic understanding, however, the systematic uncertainty in GOP estimates will remain poorly known.
Even small errors in the λ 17/18 and θ 17/18 values for respiration have outsized consequences on triple-isotope-based calculations of GOP. 22,58 As can be seen in Figure 1b This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.  22  Clumped-isotope geochemistry uses D n notation in lieu of d notation because it focuses on how isotopes are distributed among all possible isotopologues of a molecule. 62 The number of This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
12 isotopologues containing more than one rare-isotope substitution is compared to that expected for a random (stochastic) distribution of isotopes. (12) where (13) .
This key difference between these definitions and similar-looking ones above (e.g., Eq. 1) is that the ∆ 36 value is explicitly normalized against the sample's δ 18 O value. As written, then, ∆ 36 variations are independent of bulk isotopic composition and represent a new set of constraints on biogeochemical isotopic fractionation.
At this point it is useful to separate the processes that fractionate isotopologues into two broad categories: physical processes, which may move molecules but preserve their bonds, and chemical processes, which may alter, break, or make new bonds.
An example of a physical process that preserves bonds is diffusion. Knudsen diffusion, where gas flows through an orifice smaller than the mean free path of that gas, yields a calculable set of θ 33/34 , θ 35/34 and θ 36/34 values that can be compared directly to experiments 28 . The flux of gas through the orifice is proportional to its velocity, and therefore inversely proportional to the square root of its mass. The relevant fractionation factor a can thus be predicted to be (14) This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. 13 where m is the molecular mass of the isotopologues. Similar to Eq. 8 63 29 Chemical processes tend to fractionate over a different range of θ values. Here, we use an instructive example, the kinetically controlled breaking of an O-O bond. Equation 10 can be adapted to describe this process if (1) the fractionation is dominated by the atoms moving apart along the reaction coordinate and (2) isotope effects on activation energy are negligible. In that case, using the reduced masses of O 2 , μ, instead of atomic masses in eq. 10 provides an approximate description of the θ values. 33 To make this distinction clear, θ values will be denoted as θµ. Reduced masses are calculated as This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
14 where M and m are the masses of the two atoms participating in the bond (e.g., both 17.99916 amu for 18 O 18 O). Therefore, (17) .

Respiration experiments
Surface waters were sampled from Lake Houston's Deussen Park (29°55'7.1256'' N, 95°8'56.9076'' W) on June 15 and July 5, 2016 (hereafter referred to as LH1 and LH2) by casting a 20L Nalgene container from the dock and allowing it to fill with water. Upon retrieval, water was immediately partitioned into 300mL glass Wheaton bottles, stoppered and covered in This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
15 aluminum foil to suppress light-dependent reactions in favor of dark respiration. Wheaton bottles were prepared for this experiment by acid washing in 50:50 HCl and 50:50 NaOH cold rinses, then autoclaved. Evacuated 2L flasks pre-poisoned with 400 µL of a saturated HgCl 2 solution were filled with ~400 mL of cast water, resulting in a final HgCl 2 concentration in excess of that recommended to halt respiratory activity in order to determine the starting conditions (t 0 ) of the following respiration experiment. 64 Waters were transported to Rice University's Stable Isotope Lab for incubation at room temperature (measured as 25 ± 2°C); the time elapsed from first sampling at Lake Houston to lab was ~2 hours. Wheaton bottles were placed in a dark cabinet out of direct sunlight to limit temperature variations for the remainder of incubation (up to ~440 hours). One bottle was immediately unsealed and monitored with a calibrated Lazar Micro Oxygen Electrode to guide sampling throughout the experiment.
At intervals during the respiration experiment, Wheaton bottles were unsealed and water siphoned into pre-evacuated and poisoned 1L, 2L, and 5L flasks to halt respiration. The size of the flask was determined by the number of Wheaton bottles needed for combination to keep the final amount of O 2 to be analyzed between ~60-90 µmol and to keep the total volume of water less than half that of the flask. Filled flasks were transferred to an orbital shaker, where the dissolved O 2 equilibrated with the flask headspace for at least 48 hours before purification.

O 2 extraction, purification and mass spectrometry
Equilibrated flasks were inverted and the water was pumped through an evacuated flask until <1 mL water remained in the poisoned flask. At this point, the poisoned flask was isolated and moved to a vacuum line where the residual water was held at -40°C and headspace air was transferred, This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
through two U-shaped traps held at -196°C (for drying), into a silica gel finger held at -196°C for 30-45 min. Subsequent sample purification and mass spectrometry methods have been described previously.27, 65 Briefly, the samples of dissolved air on silica gel fingers are moved to an automated gas-chromatograph (GC)-based O2 purification system. Gas samples are expanded off of the silica gel in a room-temperature water bath, cryogenically dried once more, then preconcentrated onto a silica-filled U-trap prior to GC injection. Upon injection, the U-trap is heated to 90°C and flushed with helium for transfer to the GC column, which is held at -80°C to separate O2 from Ar and N2. Baseline resolution achieved by cryogenic separation of O2 and Ar permits This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. each peak to be integrated (see Figure 2). Then, the ratio of these O2 and Ar peak areas are corrected for sample-size nonlinearity and compared to those in air to obtain δO2/Ar values: (18) These values are reported as per-mil deviations from air and show long-term reproducibility of ±4‰ (1σ). After the Ar is eluted, O 2 is recollected on a U-trap filled with silica gel and held at -196°C until the collection is finished and excess He is pumped out. Finally, O 2 is warmed to 90°C for 5 minutes prior to transfer to the isotope-ratio mass spectrometer (IRMS).
After this final purification step, O 2 is transferred cryogenically to a cold prep finger adjacent to the mass spectrometer sample bellows for expansion and measurement relative to a laboratory working gas. The mass spectrometer used is a high-resolution Nu Instruments Perspective IS IRMS that typically achieves a mass resolving power (MRP) over 4000 where MRP = m 36 /(m 95%m 5% ). This increase in MRP from prior methods resolves mass interferences that were known to be an issue such as Ar in addition to others that were not. 28 For example, 35 This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
18 atmospheric O 2 dissolved in water. 28,29,63,66 In addition, bulk isotope compositions are corrected for pressure-baseline effects that yield a triple-oxygen difference between atmospheric O 2 and San Carlos Olivine that has been reproduced across three independent labs within 10 ppm. 27,30,36,67 This result implies that our measured λ 17

Computational methods
Equilibrium isotope effects for chemical reactions can be calculated using reduced partition function ratios β for chemical species according to the formula . 41,68 In the oxygen system, the formula computes the isotopic fractionation relative to atomic O vapor. It can be conceptually split into three parts indicated by the subscripts: the classical massand moments-of-inertia (MMI) term, and the quantum-mechanical zero-point energy (ZPE) and excited-state (EXC) terms. The product comprising β is across all harmonic vibrational modes i.
The quantity U i = hcυ i /k B T is the energy associated with any one vibrational mode, where h is Planck's constant, c is the speed of light, k B is Boltzmann's constant, T is the temperature, and υ i is the harmonic vibrational frequency in cm -1 . The asterisk denotes a heavy-isotope substitution.
We have omitted the symmetry numbers in Eq. 20 because each isotopic species is considered explicitly in our calculations (see below).
First-principles calculations were performed on several analogues of biological heme-enzyme active sites. Density functional theory (DFT) was used to estimate the equilibrium isotopologue This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
fractionation associated with reversible heme-O 2 binding at Cytochrome c oxidase (CcO) and other simplified oxyheme binding sites (Fig. 3). The reaction is treated as an isotope exchange reaction between free and bound O 2 , i.e., (20) The calculations thus implemented simulate the isotopic effects of reversible O 2 binding with ferrous iron, a scheme exhibited by a number of important enzymes. 69,70 The approach may be quantitative for θ values if distal electron-withdrawing moieties in enzymes have similar effects on equilibrium isotopologue effects. 71 Cytochrome c oxidase, for example, binds to O 2 in an endon configuration with distal copper instead of the distal histidine residue found in hemoglobin, but the bonding scheme may be similar. 70 This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.

20
Three levels of theory were used to study heme-O 2 binding site models: (1) the Becke 3parameter Lee-Yang-Parr functional using the D3 version of Grimme's empirical dispersion with Becke-Johnson damping (B3LYP-D3BJ), [77][78][79] (2) the one-parameter hybrid Tao-Perdew-Staroverov-Scuseria functional (TPSSh), 80,81 and (3) the second-order generalized gradient approximation functional of Pervati, Zhao, and Truhlar (SOGGA11). 82 Each level of theory is paired with the Weigand-Ahlrichs split-valence Gaussian basis sets def2-SVP and def2-TZVP, which contain polarization functions and shows good accuracy for d-block elements. 83 The first two functionals are global hybrid functionals, meaning that they treat the electron exchange-correlation energy in the Kohn-Sham DFT formulation by including a small amount of Hartree-Fock ("exact") exchange energy; B3LYP optimizes the exchange term using three parameters to yield ~20% exact exchange, while TPSSh uses one parameter to include 10% exact exchange. For metal-ligand systems, these semiempirical modifications compensate for the tendency of nonempirical functionals to bind ligands too strongly. 84 While B3LYP tends to bind ligands to d-block elements too weakly 84 , DFT studies of heme-O 2 systems suggest that empirical dispersion corrections can mitigate these shortcomings. 71,85 In contrast, TPSSh was been found to be reasonably accurate in heme-O 2 systems without the need to invoke additional empirical corrections. 71,86 The SOGGA11 functional is a pure, non-hybrid density functional based on the Perdew-Burke-Ernzerhof (PBE) functional hat has its exchange-correlation formulation optimized against a suite of molecular training sets for broad chemical applicability. 82,87,88 For isotope fractionation calculations, the accuracy of harmonic vibrational frequencies is paramount. A recent assessment of optimal harmonic frequency scaling factors suggests that both the B3LYP method yields accurate frequencies without need for scaling. 89 The TPSSh functional has not been subject to such a study, but its scaling is likely intermediate between the This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.

21
TPSS (0% exact exchange, harmonic scaling factor 1.0194 when paired with def2-TZVP) and TPSS0 (25% exact exchange, harmonic scaling factor 0.9871 when paired with def2-TZVP basis) functionals. 89 Optimal frequency scaling factors for SOGGA11 paired with the def2-TZVP basis set have not been reported. For mass-dependent fractionation, cancellation of errors may limit errors related to methodological bias regardless of the nature of the frequency error. 90 The use of three DFT methods is a test of this assertion: one expects that the range in predicted single-isotope (i.e., 18 O/ 16 O) fractionation will be larger than for multi-isotope/isotopologue trends because absolute harmonic frequency errors are cancelled during calculation of θ values.
Anharmonic corrections were not included due to computational expense; previous work suggests they contribute ~0.001 errors in θ values. 91,92 Three molecular models were investigated: a neutral free oxyheme (Heme-O  Calculations were performed in the Gaussian16 electronic structure program, revision B.01. 94 To simulate the effects of a solvent dielectric, a conductor-like polarizable continuum model (CPCM) for water (dielectric constant ε = 78) was used for Heme-O 2 and Heme(+His)-O 2 structures, while a CPCM with ε = 4 was used to simulate the CcO active site model being embedded within a supramolecular protein structure. 95 We note that while the choice of dielectric constant did not affect the results qualitatively, the inclusion of explicit waters near the binding sites could have an effect on calculated isotopic fractionations; effects on θ should be muted nevertheless. Geometries of all structures were first optimized using calculated analytical force constants, after which the frequencies calculated using the same level of theory and basis set.
The presence of a stationary point was confirmed in all cases. This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
For Heme(+His-O 2 ), the distal histidine was represented as an imidazole and needed to be fixed in place relative to the porphyrin ring. As such, two calculations were performed: beginning with a previously published structure that was based on the O 2 binding site in human hemoglobin, a set of calculations-geometry optimization and frequency-was performed in which all 20 carbons in the porphyrin ring and 6 atoms of the distal imidazole were frozen (all three C atoms, the N atom, and two H atoms furthest from O 2 ). 75 Then, a second set of calculations was performed after "unfreezing" the 8 carbon atoms closest to the O 2 moiety and freezing only the nitrogen atom and the three distant hydrogen atoms of the distal imidazole.
Because of these additional constraints, the first structure had 93 of 165 possible degrees of vibrational freedom (i.e., 56% of the total), whereas the second structure had 123 (75%). Both calculations yielded similar results, and only those for the latter, which has the least truncation of vibrational modes, are reported here. This approach resembles that employed by Rustad and coworkers to determine accurate isotopic fractionation factors from mineral clusters, wherein atoms more than two bonds away from an isotopic substitution were assumed not to contribute significantly. 96,97 Equilibrium isotopologue fractionation between the heme-O 2 adduct and free (CPCM- This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
23 included to account for the subtle isotopic preference for the heavy isotope to be bound to Fe in all cases e.g., Fe- 18   This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.

Bulk and clumped isotope variations for Lake Houston respiration
Values of δ 18 O, Δ´1 7 O, Δ 35 and Δ 36 for LH1 and LH2 are listed in Table 1  increase from 1.03‰ to 1.3‰ (Figure 5c). Values of Δ 36 increase from 1.77‰ to 1.92‰ for LH1 and from 1.95‰ to 2.58‰ for LH2 (Figure 5d). This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. Qualitatively, these results are consistent with those associated with respiration during a previous multi-phase terrarium experiment. 57 The previously observed "anomalous" covariation This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.  This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
27 regressing ln(R/R 0 ) vs. ln f to obtain α -1, where f is the proportion of initial O 2 remaining in the sample (Figure 6a) determined by assuming that dO 2 /Ar t=0 is 100% of the O 2 .
The results are robust: removing any one data point from each time trace yields best-fit 18 α resp and θ resp values within the stated uncertainty limits. Because these experiments may not have sampled identical microbial communities, some heterogeneity is expected with respect to both This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
28 initial state and effective fractionation factors. However, the large water volumes required for analysis seem to limit such effects.   This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
Bond lengths and harmonic vibrational frequencies (ν) of key bonds in the oxyhemes studied ( Fig. 3) are listed in Table 3

Some spin contamination was identified in the hybrid DFT wavefunctions (likely arising from
Hartree-Fock exchange), but previous work on oxygenated metalloprotein analogues suggests that such has little effect on the optimized geometries. 100 Bond lengths derived from the SOGGA11 functional (which lacks exact exchange) showed better agreement with experimental crystal structures, but solvated structures appear to be in better agreement with the TPSSh and B3LYP-D3BJ structures. 101 A full investigation of the accuracy of these methods is beyond the scope of this work in part because the associated errors seem to largely cancel when θ values are computed (see below). This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. cm -1 for CPCM-solvated O 2 ) showed the best agreement with experiment (e.g., 1580 cm -1 for O 2 gas 102 ), followed by TPSSh and then B3LYP-D3BJ. We note that O 2 dissolution is associated with a frequency shift of ~9 cm -1 , which is not captured in these calculations. 106 Nevertheless, the resulting fractionation is small, contributing errors of <1‰ in 18 ε and negligible errors in θ values. 29,98 The best point of comparison for 16 O- 16 O stretching frequencies in the oxyhemes (given a lack of experimental data) is between Heme(+His)-O 2 and oxyhemoglobin, which have vibrational frequencies of ~1200 cm -1 and 1155 cm -1 respectively. 103 We consider these to be in good agreement considering that the observed vibrational frequency includes anharmonic effects. Calculated isotopologue fractionations between 5°C and 50°C are shown in Fig. 7, and values at 25°C are listed in Table 4. All 18 ε values are predicted to increase in magnitude with increasing temperature, although the total change over the temperature range is less than 1‰.
This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
Results from calculations utilizing the def2-SVP basis set yielded 18   This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
34 nearly negligible variation between levels of theory and basis suggests that these results are robust for these molecular models, perhaps due to a cancellation of errors in Eq. 18. Calculated Again, SOGGA11 shows the largest deviations from other levels of theory, perhaps due to the aforementioned tendency of nonempirical functionals to overbind ligands. 84

Comparison of new laboratory experiments with previous results
Both of our incubations of natural waters yielded effective θ 17/18,resp values of 0.520, which is higher than the range of values (0.510 -0.515) reported for dark respiration. 26,44 Possible origins of these discrepancies are diverse.
One potential source of the discrepancy is pathway-specific differences in oxygen consumption. In natural waters, respiration no doubt progresses through a combination of mechanisms, even when incubated in the dark. However, an incubation of natural waters from Lake Kinneret 44 yielded a 18 ε value indistinguishable from what we observed (-21.6 ± 0.1‰ vs. -20.0 ± 1.3‰, 1σ) and a θ 17/18,resp value indistinguishable from those later obtained for pathwayspecific experiments analyzed in the same lab. 12,44 Diffusion in water has 18 ε ≈ -3.3‰ and θ 17/18 = 0.517 ± 0.002 (1σ), so diffusion control is unlikely in our experiments and cannot have driven θ 17/18,resp to higher values. 29,108 Light-dependent O 2 consumption mechanisms may have higher This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. A second potential source of the discrepancy is the difference in analytical method. Unlike other labs, we analyze pure O 2 without Ar, and our measurements are calibrated against multiple standards that together show good agreement in triple-oxygen composition relative to two other high-precision labs utilizing identical sample preparation methods. 27 Laboratories that do not separate O 2 from Ar must empirically correct for the presence of Ar using calibrated O 2 /Ar mixtures. We have shown recently that varying amounts of isobaric interferences can render triple-isotope measurements inaccurate: undesired isobars that scatter onto the O 2 analyte masses in a mass spectrometer can exacerbate "pressure baseline" scale distortion effects, which tend to decrease measured θ 17/18 values relative to their true values. Errors in θ 17/18 of 0.005 were shown to be possible using a realistic range of isobaric signals. 27 Whether these errors were present in previous experiments is not known, but our reported measurements include explicit corrections for these effects.
Regardless of the reason for the θ 17/18,resp discrepancy between labs and experiments, a discussion of theoretical context of all θ resp values (θ 17/18,resp , θ 35/34,resp , and θ 36/34,resp ) is warranted to This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
36 examine whether reported θ resp values are consistent with conceptual models for respiratory isotopologue fractionation. The novel determination of θ 36/34,resp and θ 35/34,resp values for dark respiration in this study provides additional constraints on these conceptual models and the associated θ 17/18,resp value(s) expected.

Understanding natural variability in θ values in simplified cases
In this section, we will explore how processes governed by the mass of O 2 molecules or the mass of oxygen atoms participating in O-O bonds may influence mass-dependent fractionation during respiration. Our goal is to understand the mechanistic origins of the θ resp values observed in Lake Houston waters. We will restrict our discussion to θ values because they appear to be less variable than a values for respiration, both in our experiments and in the literature. 3,12,44 Theoretical reasons for the relative invariance in θ values among different processes have been reviewed elsewhere. 32,33 Diffusive fractionation (Table 5, Rows C. D, and E) is likely relevant for respiration because O 2 must travel from the lipid membrane (where it is several times more soluble than in the surrounding medium) 110 to the respiratory enzyme. Moreover, O 2 may also need to diffuse through protein channels toward the active site (in Cytochrome c oxidase, for example). 72 If O 2containing water clusters are the relevant diffusing unit (e.g., whole clusters diffusing amongst themselves), then one can use the total masses of the isotopically substituted clusters in eq. 14, and calculate θ values via eq. 15 (Table 5, Row D). If O 2 itself diffuses through the medium rather than a larger cluster, then the process is represented using the reduced mass of O 2 and its average diffusive partner instead (Table 5, Row E). Regardless of the specific diffusive mechanism, neither of these scenarios alone are sufficient to describe the empirically determined This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
θ resp values; they are all too low, and would drive ∆′ 17 O, ∆ 35 , and ∆ 36 values down over the course of the experiment, contrary to what is observed (Fig. 5). We will examine the plausibility of this scheme with a more detailed mechanistic analysis of biological O 2 consumption below.

The role of electron-transfer processes for respiratory O 2 fractionation
Aerobic respiration is often viewed simplistically through the lens of kinetic mass-dependent fractionation of oxygen isotopes (e.g., eq. 10), but the view is conceptually flawed: O-O bondbreaking at respiratory enzymes requires the transfer of two electrons (out of four total for O 2 reduction to water), each of which weakens the O-O bond. In Cytochrome c oxidase, it occurs in two distinct steps that are shown schematically below: 72,111 (24) This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. (25) .
The first elementary step involves electron transfer from the proximal heme to O 2 ; it is reversible, so may exhibit equilibrium fractionation described by isotope-exchange reaction 21. 112 The second step is coupled with proton transfer. 74,111 Previous studies suggest it may be irreversible in many respiratory oxidoreductases, serving to "commit" the enzyme to O 2 reduction. 113  This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
Small deviations in these vibrational frequencies (up to tens of cm -1 ) do not alter the basic observations for θ values reported here.
The second electron transfer step, reaction 26, is not well understood. It is likely kinetically controlled, but the relevant transition-state structure is poorly known, preventing a confident DFT-based estimate of its kinetic isotopologue fractionation to be made at present. 74,111 Nevertheless, one can examine proton-transfer and electron-transfer limitation qualitatively. For the former case, formation of an O-H + bond is limiting; because the reaction coordinate involves bonds other than the O-O bond, the mass dependence of the reaction may lie somewhere between θ kin and θ μ values in Rows A and F of

Two-step model for isotopic fractionation at a respiratory oxidase
To quantify how these chemical-and diffusion-dominated models differ, and to determine whether the clumped-isotope data can distinguish between them, we utilized the two-step This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.

41
The first step represents O 2 diffusing to the enzyme active site. The second step encompasses both the O 2 binding and reduction steps. We chose to combine the two latter steps for this exercise because they show a similar range of θ values. Our goal is to use the model to distinguish between diffusive and electron-transfer limitation of respiratory O 2 consumption. For this reaction sequence, the measured isotopic fractionation factors (α measured ) can be described using fractionation factors for aqueous diffusion and enzymatic reduction (α diffusion and α enzyme , respectively) and a reversibility parameter r ranging from 0 to 1, where r = 1 indicates complete reversibility in the first step, i.e., no diffusive fractionation: (27) For diffusive fractionation, we use 18 Fig. 8.
This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review. We find that our experimental data are generally consistent with natural respiration expressing some diffusion limitation (i.e., implied r > 0.4). Nearly all the isotopologue fractionation factors can be explained simultaneously if the combined 18 ε enzyme value is between -30‰ and -40‰: Fig.   8 shows all our measurements are within 1s of this envelope with one exception. The exception is θ 35/34,resp in experiment LH1; in that experiment, that the expected change in ∆ 35 values is near the limits of detection, so the uncertainty in the experimental θ 35/34,resp value may be underestimated. Moreover, a kinetic 18 ε enzyme value of this magnitude is plausible as the cumulative fractionation for binding and reduction, or for an outer-sphere quantum-mechanical electron-transfer step to O 2. 114 The specific mechanism that applies to this system, however, awaits further investigation. Interestingly, the covariation of 18 ε measured and θ 17/18,resp with changes in temperature reported by Stolper et al. (2018) can be explained qualitatively by this two-step model if the diffusion limitation becomes less important at higher temperature, where molecular transport is expected to be faster. This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.
In contrast, the lower range of θ 17/18,resp values can only be explained by the model if the reaction occurs closer to the diffusion limit, e.g., r ~ 0.2 or less, and the intrinsic enzymatic 18 O/ 16 O fractionation is much stronger, e.g., 18 ε enzyme ~ -70‰. To our knowledge, such large 18 O/ 16 O fractionation has not been observed before in a non-photochemical system. Well-stirred experiments with isolated enzymes may shed light on which of these cases is more accurate.  18,55,118 However, the low-O 2 endmembers had never been observed, rendering this interpretation somewhat speculative. Previous measurements were made on O 2 + Ar mixtures, which may have been susceptible to "pressure baseline" errors and variability therein, so they may need revision. This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.

Conclusions and implications for determinations of primary productivity
Our results raise questions about the accuracy of the "canonical" θ 17/18,resp range near 0.515. The Lake Houston incubations, at face value, suggest that the natural θ 17/18,resp value could be 0.004 higher. If true, they imply that existing GOP estimates could be ~40% too high (see Fig. 1).
Downward revisions in GOP of this magnitude may resolve lingering disagreements between the triple-oxygen method and H 2 18 O incubations and disagreements between O 2 production-to-C fixation ratios observed in situ and in cellular physiology studies. 19, 21, 47, 119 They may also change estimates of ecosystem carbon export efficiencies, which are especially important in oligotrophic regions. 54,120 However, there may be fortuitous cancellation of errors in the calculation of GOP because the triple-oxygen composition of seawater typically used in those calculations may also be erroneous. 67 Indeed, revising the photosynthetic O 2 endmember to reflect more recent determinations of the triple-oxygen composition of the Vienna Standard Mean Ocean Water-2 (VSMOW2) standard could mitigate the impact of a +0.004 revision in θ 17/18,resp on GOP estimates almost entirely (Fig. 9). 30,31 To what extent marine and global GOP estimates need to be revised, then, will not be known until accurate θ 17/18,resp values, and how they change under different environmental conditions, can be obtained. The reported sensitivity of θ 17/18,resp values to temperature is especially important in this context. 26 Moreover, the use of triple-oxygen isotope variations in various archives to assess past biosphere productivity and to infer the presence of biological activity relies on accurate θ 17/18,resp values and a strong understanding of the factors that modulate its variability. 23-25, 122, 123 Thus, our work highlights the need to re-evaluate both θ 17/18,resp and the triple-isotope This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.  This is a non-peer reviewed pre-print submitted to Earth ArXiv. This manuscript has been submitted to ACS Earth and Space Chemistry for peer review.