A modified method to calculate dual-isotope slopes for the natural attenuation of organic pollutants in the environment

Two-dimensional compound-specific isotope analysis has become a powerful tool for distinguishing reaction mechanisms. Lambda (Λ), an essential and important parameter for processing two-dimensional isotope fractionation data, exhibits values specific to a reaction mechanism. In the present article, we modified the existing algorithms for calculation of lambdas based on a review of current methods. Specifically, by regressing [(1000+δE0,2)*(n1*x2)*ΔδEbulk,1] versus [(1000+δE0,1)*(n2*x1)*ΔδEbulk,2] by the York method, a novel method was developed to calculate Λs. The improved method eliminates both the influence of the nonreacting position and the initial isotope signatures. Furthermore, this method retains the advantages of a two-dimensional isotope plot, which eliminates contributions from commitment to catalysis, does not require determination of the fraction of remaining substrate, and can be constructed even from field data. Additionally, the one-sample t test is applied to generate a 95% confidence interval of the dataset of Λris for various reaction mechanisms. The ranges of 5.67–24.8, 8.54–9.80, 0.51–8.35, 25.2–36.8, and 7.09–21.9 are applicable for the oxidation of C-H bonds (ZC=1, ZH=3; ZC and ZH are the number of indistinguishable carbon and hydrogen atoms in intramolecular competition, respectively), oxidation of C-H bonds (ZC=1, ZH=4), aerobic biodegradation of benzene (ZC=6, ZH=6), methanogenic or sulfate-reducing biodegradation of benzene (ZC=6, ZH=6), and nitrate-reducing biodegradation of benzene (ZC=6, ZH=6). The accumulation and correction of these values will make the data measured in the field easier to interpret.


Introduction
Compound-specific isotope analysis (CSIA) has been extensively explored as a useful tool to identify contaminant sources and monitor the extent of pollutant degradation (Elsner 2010;Elsner and Imfeld 2016). By measuring the kinetic isotopic effects (KIEs) for compounds undergoing biochemical reactions, CSIA can also provide insight into biodegradation reaction mechanisms (Elsner et al. 2007). KIEs for a specific reaction mechanism are fully expressed only if the dominant isotope fractionation step is the slowest step, known as the rate-determining step or rate-limiting step (Jaekel et al. 2014). However, the isotope-sensitive step (i.e., the bond conversion) is occasionally preceded by a non-or slightly fractionating process, such as transport to and adsorption at reactive sites, or the formation of enzyme-substrate complexes, which will make the measured apparent kinetic isotope effect smaller than the intrinsic KIE Feisthauer et al. 2011). The masking factors often include high microbial cell densities, low substrate or electron acceptor bioavailability, substrate transport across cellular membranes, and the reversibility of enzymatic reactions (Jaekel et al. 2014). The expression "commitment to catalysis" (Northrop 1981) was introduced to describe this influence. To mitigate the masking effects of biochemical reactions, the simultaneous measurement of isotope signatures of two elements involved in isotope-sensitive steps (two-dimensional CSIA) has been proposed (Fischer et al. 2007;Jaekel et al. 2014). The slope of a dual-isotope plot, lambda (Λ), can be specific to one certain reaction mechanism.
Generally, four methods are used to calculate Λs, as follows: (i) Calculating Λs from the fractionation factor at the reacting position (α rp ) . When there are heavy isotopes in the reaction site of the molecule, α rp is usually masked by the commitment to catalysis (C) and the influence of C can be described by Eq. ((1) , where the superscripts "⊂" and "⊄" represent the presence or absence of commitment to catalysis. As masking effects always equally affect the two elements involved in biochemical reactions Zwank et al. 2005), the combination of two-dimensional isotope equations elegantly eliminates contributions from commitment to catalysis, as shown in Eq. ((2).
(ii) Calculating Λs and its uncertainty by dual-isotope plots that determine regression with ordinary linear regression (OLR). The rate-limiting step of nonfractionation, including the case in which the substrate is absorbed by cells or combined with an enzyme (caused by the enzyme masking effect), has an equal effect on the isotopic effects of different elements. Therefore, the masking effect of commitment to catalysis can be reduced to a certain extent with the ratio calculation. Isotope signatures of different elements (δE 1, t , δE 2, t ) can be plotted on the horizontal and vertical coordinates respectively and the slope of ΔδE 1, t versus E 2, t (ΔδE t = δE t − δE 0 ), denoted by Λ OLR , is obtained from OLR.
(iii) Calculating Λs and their uncertainties by dual-isotope plots that use regression with the York method (Ojeda et al. 2019). In plotting dual-isotopes, OLR only minimizes the sum of squares of the isotope signature set as the y-variable, without considering the measurement error of the other variable. Measurement errors are inherent to both variables and neither of them can be ignored. Λ york is obtained by regressing measured data through the York method, which incorporates uncertainty in both isotope measurements (Ojeda et al. 2019); Λ york was proposed as a more appropriate estimation of lambda than Λ OLR .
(iv) In addition to the above, there is another way to calculate Λ (Eq. ( (3)) (Jaekel et al. 2014): where ε rp,1 and ε rp,2 are isotopic enrichment factors for the two isotopes in the reactive positions.
To the best of our knowledge, there is limited comprehensive and in-depth comparison among these four methods for data processing. Nevertheless, there remain puzzling disparities among the Λs calculated with the four methods. The objective of this paper is to propose a modified method for calculation of Λs that can easily be used to interpret field survey data.

Data collection
Published studies using dual isotope (C, H) fractionation allowing practitioners to identify transformation mechanisms are far more common than studies using other isotopes. This provided a good opportunity to evaluate the uncertainties associated with the current mathematical m o d e l s w i t h c a r b o n a n d h y d r o g e n i s o t o p e s . Experimental data were collected from previous studies involving biodegradation of organic contaminants, including carbon and hydrogen isotope measurements (Bennett et al. 2018;Cui et al. 2017;Elsner et al. 2007;Feisthauer et al. 2011;Fischer et al. 2008;Holler et al. 2009;Hunkeler et al. 2001;Jaekel et al. 2014;Kinnaman et al. 2007;Mancini et al. 2008;2003;McKelvie et al. 2009;Rasigraf et al. 2012;Vogt et al. 2008;Zhang et al. 2019).

Data analysis
All calculations of Λs were performed in R using the native functions for OLR and the IsoPlotR package for York regression, and the pseudocode for the regressions was provided in the prior literature (Ojeda et al. 2019). Sensitivity analysis using Crystal Ball 11 software was conducted to reveal the key input factors that primarily affect the estimation of the new Λ. The contribution of each input variable to the output (the value of Λ) variance was assessed, and the sensitivities of each input variable relative to the others were examined.

Method improvement
The merits and drawbacks of existing methods for calculation of lambda were reviewed carefully. Then, the wheat was winnowed from the mountains of chaff, and an improved method for lambda calculation was developed, mainly by working hard to eradicate the influences of the initial isotope signatures and nonreacting positions based on those current methods. Using the modified algorithm (Λ ri ) to elucidate the reaction mechanism, the data collected in the prior literature were recalculated for determination of Λ ri , and one-sample t test was used to generate the 95% confidence interval of the range of Λ ri for certain reaction mechanisms.

Evaluation effects
By weighing the relative merits of the modified method first and comparing the lambda values calculated with different methods, the validity of the new mathematical model was determined. Specifically, five algorithms, including four existing algorithms and one derived here, were used to calculate lambda. All the calculated results were used to evaluate the five methods, and especially the new approach (Table 1 and Fig. 1).

Merits/demerits of algorithms for lambdas
Mechanisms of various biochemical reactions can be compared in a straightforward manner because distinct mechanisms usually lead to diverse Λ α s. In addition, the major advantage of this approach is the elimination of the masking effect caused by commitment to catalysis. However, this method is rarely used in the literature (Rasigraf et al. 2012). First, the error in Λ α s may increase after many operations with the original data. An example can be found in a previous study (Feisthauer et al. 2011), in which the uncertainty ranged from 1.2 to 3.3, with 11-33% uncertainty in the calculated value of Λ α (8.2-10.5). Furthermore, determining the precise fraction of the remaining substrate ( f ) is necessary to calculate α rp through Eq. ((2). Actually, measurement of f is very difficult in field investigations. Moreover, this method may be inappropriate for isotopes with high natural abundance (e.g., S and Cl) because it neglects the effect of molecules with more than one heavy isotope for isotope fractionation with a limited number of molecules (or low abundance molecules).

Lambda OLR (Λ OLR )
The major advantage of Λ OLR is that a two-dimensional isotope plot can easily be constructed, even from field data, without measuring the fraction of remaining substrate. Second, any two isotopes (even Br (Tang and Tan 2018) or Cl (Rodriguez-Fernandez et al. 2018)) can be used to obtain Λ OLR without the constraint of natural abundance. Third, the determination can be completed in a one-step regression, thus reducing the transmission of errors. Last but not least, because Λ OLR is approximately equal to Λ α (Elsner et al. 2007;Vavilin and Rytov 2016) (the derivation is in the literature (Elsner et al. 2007) and revised below), distinct mechanisms usually lead to different Λ OLR values (Braeckevelt et al. 2012). Considering the above advantages, Λ OLR has been widely used in cases involving two-dimensional isotopes (Badin et al. 2014;Solano et al. 2018) (or three (Kuder et al. 2013, Van Breukelen et al. 2017).
However, there are still some flaws in this method. First, according to Eq. ((1), the masking effect of commitment to catalysis on isotope fractionation is not linear. Therefore, it could not be eliminated exactly by simply using the ratio calculation. Especially in the case of high isotope fractionation (e.g., ε>100‰), which is often observed with hydrogen isotopes, the calculation of Λ OLR is less reliable than that of Λ α (Feisthauer et al. 2011). Furthermore, calculating Λ OLR by OLR of ΔδE 1 versus ΔδE 2 has not been corrected with respect to nonreactive locations. This occasionally makes Λ OLR deviate from Λ α considerably. Toluene degradations by several anaerobic cultures were investigated ) by twodimensional CSIA. Λ OLR and Λ α values are shown in Table 1 and Fig. 1. Clearly, Λ OLR is more than twice as large as Λ α , which means that Λ OLR values cannot be compared with those produced with other molecules in the same reaction type. However, there are only a few restrictions for Λ α , and one could be compared with those for other molecules reacting via a similar reaction mechanism . As a result, the relationship between ΔδE and Λ α is reduced and Λ ri and other factors that may affect the accuracy of the dual-isotope plot are considered.

Lambda york (Λ york )
Compared with Λ OLR , Λ york has some advantages, such as regression symmetry (Ojeda et al. 2019). Similar to Λ OLR , nevertheless, the deviations between Λ york and Λ α are pronounced in some cases (Table 1 and Fig. 1), which indicates that Λ york needs to be further considered.

Lambda rp (Λ rp )
Λ rp takes the influence of nonreacting positions into account and is much closer to Λ α than to Λ OLR or Λ york (Table 1 and Fig. 1). However, this equation does not make much sense. On the one hand, Λ rp does not completely eliminate the masking effect of commitment to catalysis, unlike Λ α . The relationship between Λ rp and Λ α is shown in Eq. ((4) (the derivation is provided in the Supplementary information).
On the other hand, Λ rp has almost the same drawback as Λ α , in that the precise fraction of the remaining substrate must be determined and more calculation steps lead to greater uncertainties. Therefore, it is wiser to calculate Λ α by Eq. ((4) than by ε rp, 1 /ε rp, 2 , if ε rp, 1 and ε rp, 2 are known.
In conclusion, Λ α and Λ rp are more suitable for laboratory data and mechanism-oriented analysis because f is easy to measure, the influence of nonreacting positions is eliminated, and different compounds reacting via the same mechanism can be compared. Λ OLR and Λ york are easier to calculate from field data and usually differentiate mechanisms for only one specific compound. Therefore, a modified method to process field data and compare reaction mechanisms for different compounds is presented below.

Results and discussion
Modified algorithm for lambda E q u a t i o n ( ( 5 ) ( d e r i v a t i o n i s p r o v i d e d i n t h e Supplementary information) is proposed to calculate lambda (it is defined as Λ ri in the present study), and the disadvantages mentioned above have been swept away. Here, the York method is selected to regress [(1000+δE 0 , 2 )*(n 1 *x 2 )*ΔδE b u l k , 1 ] versus [(1000+ δE 0,1 )*(n 2 *x 1 )*ΔδE bulk,2 ] since regression symmetry and the measurement error are accounted for in both variables, Furthermore, Eq. ((5) retains the advantages of a twodimensional isotope plot; it has less uncertainty, there is no need to determine the fraction of remaining substrate, and the plot can be constructed even from field data. To verify the new algorithm, all five algorithms were used to recalculate Λs, and the results are listed in Table 1. Λ α is set as the criterion for two-dimensional CSIA, as it has proven to be mechanism-specific . The deviation between Λ ri calculated by Eq. ((5) and Λ α is very small, even in the case of a large deviation between Λ α and Λ OLR or Λ york (Table 1 and Fig. 1). As a result, Λ ri calculated by Eq. ((5) is strongly recommended to describe the relationship between two isotopes. Some factors that may affect the accuracy of Λ ri are discussed below.

Merits of the new method
This method looks simultaneously at many aspects of uncertainty in the algorithm for lambda. Specifically, the improved method eliminates both the influences of the nonreacting position and the initial isotope signatures. At the same time, this method retains the advantages of a two-dimensional isotope plot. For example, it eliminates contributions from commitment to catalysis but does not require determination of the fraction of remaining substrate and can be constructed even from field data.
The isotope signature measured with gas chromatography/ isotope ratio mass spectrometry can only provide average bulk isotope ratios, in which the changes will be smaller than those at the reacting positions. In several molecules, there is little difference between Λ ri and Λ OLR or Λ york . In the enrichment of 13 C and 2 H during monooxygenase-mediated biodegradation of 1,4dioxane (Bennett et al. 2018), 1,4-dioxane is the target molecule, C and H are the elements considered, and n C =8, n H =8, x C =4, x H =4, so n 1 *x 2 n 2 *x 1 ¼ 1. Nevertheless, in most cases, n 1 *x 2 n 2 *x 1 is not equal to 1, and some results even deviate greatly. For toluene, the deviation is 38.1% (Table 1 and Fig. 1).
One may say that the determination of n and x will complicate the analysis. However, this step is necessary. Λs obtained from average bulk isotope ratios (Λ OLR or Λ york ) imply a hypothesis that all atoms of the considered element in the molecule are in reactive positions. On the other hand, distinct molecules usually present diverse values of n 1 *x 2 n 2 *x 1 , which greatly restricts the comparison between them even if a similar reaction mechanism is occurring. The transformation of Λ OLR or Λ york into Λ ri makes it convenient to compare distinct molecules with each other, irrespective of their molecular size. For the same biochemical reaction, Λ ri could be compared in different contaminants, even for structurally very dissimilar contaminants (Table 2 and Fig. 2).  Table 1). "Anaerobic biodegradation of toluene," "aerobic biodegradation of DMP," "aerobic biodegradation of DEP," "anaerobic degradation of propane," and "aerobic degradation of 1,4-dioxane" are abbreviated as "anaerobic toluene," "aerobic DMP," "aerobic DEP," "anaerobic propane," and "aerobic 1,4-dioxane," respectively  Table 2) In one-dimensional CSIA data processing, practitioners often ignore the effect of the initial isotope signatures but focus more on the range of measured delta values. However, different initial isotope values may lead to different conclusions. The theoretical derivation and an example are given in the Supplementary information. Similarly, the initial isotope signatures have a certain influence on the two-dimensional CSIA data processing. In general, the initial isotope signature ratios ( 1000þδE 0;2 1000þδE 0;1 ) are very close to unity and have minimal impact on Λs. The mathematical error is less than 5% when the deviation between δE 0, 1 and δE 0, 2 is less than 50‰ (|δE 0, 1 − δE 0, 2 | < 50‰). The effect is ignored when dealing with the data in the overwhelming majority of papers, and the information about initial isotope signatures provided in some papers is shown in Table S1. However, a negative situation may emerge in some cases. The changes in the carbon and hydrogen isotope signatures of ethane during aerobic biodegradation were investigated (Kinnaman et al. 2007), and δ 2 H 0 and δ 13 C 0 were approximately −32.1‰ and −347‰ respectively. In addition, the arithmetic result of (1000 + δE 0, C )/(1000 + δE 0, H ) was 148%. In other words, the impact of the initial isotope signatures may cause Λs to deviate by 48%. In addition, due to the different original sources and production procedures, the initial isotope signatures used in each laboratory and between organics vary, which limits the comparison of Λs between different laboratories and molecules. Additionally, (1000 + δE 0, 2 )/(1000 + δE 0, 1 ) is easy to obtain, and use of the initial isotope signatures is strongly recommended.

Sensitivity analysis of the new algorithm (Λ ri )
A quantitative sensitivity analysis to evaluate the variability and uncertainty of parameters in Λ ri estimation was conducted with Monte Carlo simulations (more information is provided in the Supplementary information). The sensitivity analysis results are shown in the form of percentile contributions to variance (Fig. 3). The two most influential variables that contributed to the total variance of Λ ri are ΔδH t /ΔδC t (the changes of isotope signatures) and (n H *x C )/(x H *n C ) (the influence of nonreacting positions). Although the contributions for (1000+δC 0 )/(1000+δH 0 ) are quite limited (1.9%), the influence of initial isotope signatures should not be ignored and the reason is discussed above.

Implications for reactions mechanisms
Different Λs to illustrate the attenuation mechanisms The fraction of the remaining substrate, a necessary parameter for calculating Λ α and Λ rp , is difficult to determine in field investigations, which limits the use of Λ α and Λ rp . As with recent recommendations, there has not been sufficient time for Λ york to be routinely adopted (Ojeda et al. 2020), and the difference between Λ york and Λ OLR is not significant, as shown in the results of previous literature (Ojeda et al. 2019). Therefore, the values of Λ OLR and Λ ri , two of the five different methods used to calculate Λ, are usually discussed for known mechanisms. Acid hydrolysis of fuel oxygenates (including methyl tert-butyl ether, MTBE; ethyl tert-butyl ether, ETBE; tert-amyl methyl ether, TAME) is expected to occur via an S N 1 mechanism (Elsner et al. 2007;Rosell et al. 2012) and secondary hydrogen isotope effects may be anticipated. Some values of Λ OLR much larger than the usual secondary hydrogen isotope effects were measured in the reaction (Table 2 and Fig. 2), which are statistically identical to some Λ OLR values measured for aerobic C-H bond cleavage (Feisthauer et al. 2011). In addition, some Λ OLR values determined in the aerobic biodegradation of toluene (Mancini et al. 2006;Vogt et al. 2008) were much larger than those for methane (Feisthauer et al. 2011), even though both of them undergo aerobic C-H bond cleavage. However, when Λ ri is used for data analysis, (i) Λ ri s of 1,2-dichloroethane biodegradation and fuel oxygenate acid hydrolysis are statistically identical; (ii) Λ ri s of fuel oxygenate acid hydrolysis and methane aerobic transformation are significantly distinguished and (iii) the difference in Λ ri s for methane and toluene aerobic biodegradation is decreased; there are some differences in Λ ri values resulting from the same reaction mechanism, which is mainly

Elucidation of reaction mechanisms
The relationship between Λ and KIE (Eq. ( (6)) is derived in prior literature . KIE is mechanism-specific (Ojeda et al. 2019), so Λ ri can be compared for the identical mechanism with the same Z.
where Z is the number of indistinguishable atoms in intramolecular competition. When a set of Λ ri s from certain known reactions are collected, the 95% confidence interval will be generated through a one-sample t test. If the value of an individual Λ ri measured from any unknown reaction falls within the confidence interval, the substrate is likely to have reacted with the same mechanism. Therefore, it is necessary to collect sets of Λ ri s from certain known reactions (Table S2 in the Supplementary information) and to generate the 95% confidence intervals. The ranges of 5. 67-24.8, 8.54-9.80, 0.51-8.35, 25.2-36.8, and 7.09-21.9 are applicable for the oxidation of C-H bonds (Z C =1, Z H =3), oxidation of C-H bonds (Z C =1,Z H =4), aerobic biodegradation of benzene (Z C =6,Z H =6), methanogenic or sulfate-reducing biodegradation of benzene (Z C =6,Z H =6), and nitrate-reducing biodegradation of benzene (Z C =6,Z H =6). In the case of C-H bond oxidation, some values of Λ ri are not in the range 5.67-24.8 (such as Λ ri =40.3 (McKelvie et al. 2009) and 4.19 ). Larger values may be due to the hydrogen secondary isotope effect, which causes strong hydrogen isotope fractionation (Elsner et al. 2007) and therefore larger Λ ri s.
From another point of view, the one-sample t test focuses on the overall data, not on particularly high or low values of Λ ri s that may theoretically be outliers of the data group. Nevertheless, due to the limited sample data, the range set for various reaction mechanisms may be slightly inaccurate. In future research, more data from various reaction mechanisms could be collected to generate the 95% confidence intervals of all mechanisms.

Conclusions
Sensitivity analysis of the new algorithm (Λ ri ) shows that the contributions for ΔδH t /ΔδC t , (n H *x C )/(x H *n C ), and (1000+ δC 0 )/(1000+δH 0 ) accounted for approximately 0.73, 0.25, and 0.02, respectively, of the variability and uncertainty of calculating Λ ri . Λ ri takes the influence of the nonreacting position and the initial isotope signatures into account, leading to an appropriate estimation of lambda. The ranges of Λ ri values responsible for several specific reaction mechanisms were generated using Student's t test. These ranges facilitate the assignment of the reaction mechanism in the field. Λ ri was calculated by Eq. ((5) with isotope signatures (δ t ) and assumed values of n and x. If Λ ri measured from an unknown reaction falls within the range of a certain reaction mechanism, the unknown reaction most likely proceeded via the same mechanism. This is of great significance for ecological risk assessment and the design of remediation strategies for pollutants. By increasing the number of effective substances that could promote the reaction, for example, the problem of contaminants in the environment could be greatly reduced. In conclusion, this set of best practices will help practitioners select the appropriate method to process the data of two-dimensional CSIA and apply it to the natural attenuation of organic pollutants in the environment.
Availability of data and materials All data generated or analyzed during this study are included in this published article and its supplementary information files.
Author contribution Jin-Ru Feng: investigation, data curation, writing original draft. Hong-Gang Ni: validation, supervision, writingreview and editing, funding acquisition. Both authors read and approved the final manuscript.
Funding This study was financially supported by the Shenzhen Science a n d T e c h n o l o g y R e s e a r c h a n d D e v e l o p m e n t F u n d (JCYJ20180302150410243) and the National Natural Science Foundation of China (No. 21876004).

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication Not applicable.
Competing interests The authors declare no competing interests.