Modeling Bioavailability Limitations of Atrazine Degradation in Soils

Pesticide persistence in soils is a widespread environmental concern in agro-ecosystems. One particularly persistent pesticide is atrazine, which continues to be found in soils and groundwater in the EU despite having been banned since 2004. A range of physical and biological barriers, such as sorption and mass-transfer into bacterial cells, might limit atrazine degradation in soils. These effects have been observed in experiments and models working with simplified systems. We build on that work by developing a biogeochemical model of the degradation process. We extended existing engineered system models by including refined representations of mass-transfer processes across the cell membrane as well as thermodynamic growth constraints. We estimated model parameters by calibration with data on atrazine degradation, metabolite (hydroxyatrazine) formation, biomass, and isotope fractionation from a set of controlled retentostat/chemostat experiments. We then produced site-specific model predictions for arable topsoil and compared them with field observations of residual atrazine concentrations. We found that the model overestimated long-term atrazine biodegradation in soils, indicating that this process is likely not limited by bioavailability or energetic constraints of microbial growth. However, sorption-limited bioavailability, could explain the long-term fate and persistence of the main degradation metabolite hydroxyatrazine. Future studies should seek alternative controls that drive the observed atrazine persistence in soil. This work helps to bridge the gap between engineered and natural systems, allowing us to use laboratory setups to gain insight into real environmental systems.


INTRODUCTION
The worldwide intensification of agriculture is closely linked to increased use of pesticides (Roser, 2019). Persistent pesticides are defined as those that remain in soils "in significant concentrations until the next growing season" (Craven and Hoy, 2005). Field monitoring campaigns have demonstrated the presence of residual pesticides across Europe (Silva et al., 2019).
Pesticide degradation in natural environments may be impeded by a range of physical and biological constraints. For instance, sorption of pesticides onto soil particles limits microbial access to pesticides, retarding degradation (Novak et al., 1995;Guo et al., 2000;Siek and Paszko, 2019). Moreover, spatial heterogeneity and separation of microorganisms and pesticides in soil reduces biodegradation rates (Dechesne et al., 2010;Pinheiro et al., 2018). Diffusion-limited transport across the cell membrane has been identified as a potential limiting step of pesticide degradation under low concentrations, based on observations made in engineered (chemostat and retentostat) systems (McKelvie et al., 2007;Thullner et al., 2008;Braeckevelt et al., 2012;Ehrl et al., 2018b;Kundu et al., 2019Kundu et al., , 2020. Likewise, under specific conditions, the energy produced from catabolism of some pesticides may be insufficient to support cellular energy needs, leading to pesticide persistence despite microbial accessibility (LaRowe and Van Cappellen, 2011). To date, the effect of these barriers has only been explored in the lab under controlled conditions (Ehrl et al., 2018bKundu et al., 2019) or in simulation studies based on simplified systems (Guo et al., 2000;Gharasoo et al., 2019;Marozava et al., 2019;Kundu et al., 2020).
In this work, we apply biogeochemical modeling to investigate potential factors of long-term pesticide persistence in soils. We extended existing chemostat/retentostat models  by the 1) introduction of thermodynamic growth constraints (Desmond-Le Quéméner and Bouchez, 2014;Ugalde-Salas et al., 2020) [the alternative model formulation uses a simple Monod kinetics growth ], 2) a refined formulation of mass-transfer processes across cell membranes, and 3) calibration against isotope fractionation data. These engineered systems make it possible to study the biodegradation of pesticides at growth-limiting substrate concentrations (Kundu et al., 2020), which is not possible in situ for soils. Chemostats are continuous flow systems that allow for control of steady-state microbial growth rates by adjusting substrate input and bacterial washout rates (Kuenen, 2019). Retentostats are similar devices, but retain the bacterial biomass, allowing for analysis of microbial substrate turnover at extremely low substrate concentrations (Kundu et al., 2020). We then extended the model by including equilibrium sorption and leaching in soils, and ran site-specific predictions of pesticide degradation in soil over 30 years. We compare our model predictions with residual atrazine concentration of topsoils at two study sites (arable soil) in Germany at which no atrazine has been applied for over 30 years. Albeit the long-term predictions show considerable discrepancies with the field data, our analysis provides insight into the relative contributions of model features toward long-term atrazine persistence in soils.

Model Description
Our model ( Figure 1) describes a single bacterial population (C B ) that uses atrazine (AT) as its sole carbon (C) and energy source. The core model (green background), describes behavior in engineered systems (chemostat/retentostat); it incorporates intracellular and extracellular compartments, each of which contain concentrations of both AT and hydroxyatrazine (HY). Hydroxyatrazine is produced by dechlorination of the side chain of AT. This is the first metabolic step of AT degradation. We extended the model to soil (blue background) by incorporating equilibrium sorption and leaching for each component in the extracellular compartment.

Atrazine and Hydroxyatrazine Degradation
The model describes pools of atrazine (AT) [µg L −1 ] and hydroxyatrazine (HY) in the intracellular and extracellular compartments: AT i /HY i and AT e /HY e , respectively. To take advantage of available data on isotope fractionation of AT, we split the AT pools into light (AT l ) and heavy (AT h ) isotopologues ( 12 C/ 13 C) in each compartment.
We modeled degradation of both isotopologues of AT with Michaelis-Menten kinetics, allowing for competition for binding sites. For the light isotopologue: The slightly slower degradation of the heavy isotopologue is captured by scaling the maximal degradation rate by β)1 as follows: where again k HY [d −1 ] is the maximum degradation rate, but now K HY M [µg L −1 ] is a reference concentration for growth. These two variants [Monod (M), Thermodynamic (T)] show similar behavior at high HY concentrations (such as in chemostat/retentostat systems), but differ considerably at low HY concentrations (such as in soil).

Mass-Transfer
We account for diffusive transport of AT and HY across the cell membrane where l indicates the light isotopologue, and h the heavy isotopologue, and r e [L d −1 μg −1 ] is the mass-transfer rate coefficient assumed to be the same for both compounds.

Maintenance
We incorporate metabolic maintenance requirements following the Pirt model (Pirt, 1982;Gharasoo et al., 2019): is the bacterial biomass, and Y [-] is the yield coefficient.

Input and Washout of AT, HY, Biomass
For engineered systems (chemostat/rententostat), we include a constant input of AT as: where r D [d −1 ] is the dilution rate. Additionally, we define washout terms for biomass, and AT and HY: where α [-] is 1 for a chemostat (from which biomass is washed out) and 0 for a retentostat system (where biomass is retained). The core model is described by the following system of ordinary differential equations (ODE): dAT l e dt r AT l input − r AT l e washout − r AT l mass−transfer · C B (18) FIGURE 1 | Model structure for engineered (chemostat/retentostat) systems (green) and extension for soil (blue). The model explicitly accounts for light "l" and heavy "h" isotopologues ( 12 C/ 13 C) of AT due to enzymatic transformation in the intracellular "i" and extracellular "e" compartments, as well as in the equilibrium sorption in the soil "e,S".
where f cell [µg cell −1 ] is a conversion factor from cells to carbon, and V u [L] is the volume of a single bacterium, set to 1 · 10 -15  (full details in the Supplementary Section 2). The last terms in Eqs. (16), (17) and (20) account for changes in inner cell concentrations as the total bacterial volume changes due to growth and decay.

Extension for Soil
As shown in Figure 1, we extend the core model by including equilibrium sorption and transport. We partition the extracellular concentrations of both AT isotopologues, as well as HY, into solution phase and sorbed phase concentrations: where C T [µg L −1 ] is the total extracellular concentration of AT, and HY, C L [µg L −1 ] is the solution phase concentration (AT l e , AT h e , HY e ), C S [µg kg −1 ] is the sorbed phase concentration (AT l e,S , AT h e,S , HY e,S ), θ [-] is the soil water content, and ρ [kg L −1 ] is the soil bulk density.
We relate C L and C S by the Freundlich isotherm: implemented in the model via the retardation factor: where K F (K AT and K HY for AT and HY respectively) [μg (1−nF) Kg −1 L nF ] is the Freundlich coefficient and n F (n AT and n HY for AT and HY respectively) [-] is the Freundlich exponent. Additionally, Arthrobacter aurescens TC1 and other AT degraders utilize other organic substances as C and energy source. We, therefore, assume that a minimum AT degrader biomass is maintained in soil (Klier et al., 2008): where M [µg L −1 ] is the minimum bacterial biomass in soil. Transport is restricted to convective flow: r HYe where v v [d −1 ] is the water flow per soil volume in the plough layer. We did not include abiotic degradation of AT (Kiss et al., 2007;Kiss and Virág, 2009;López-Muñoz et al., 2011), which has been observed to have a relatively small contribution compared to biotic degradation (Braeckevelt et al., 2012).
The full model for soil is described by the following system of ODEs.

Engineered Systems: Experimental Details
We calibrated two model variants (M: employing Monod-kinetics for HY degradation; T: employing thermodynamic HY biodegradation constraints) against published data from chemostat and retentostat experiments (with two replicates per experiment). Atrazine was provided as the sole C and energy source for the bacterial strain Arthrobacter aurescens TC1 (Ehrl et al., 2018b;Kundu et al., 2019). Both engineered systems were fed with an AT solution (30 mg L −1 ), with dilution rates, for the chemostat, of 0.023, 0.032, 0.048, 0.056, 0.068 d −1 , and, for the retentostat, of 0.02 d −1 . For each system at each dilution rate, concentrations of AT [µg L −1 ], HY [µg L −1 ], and living biomass [cell L −1 ] were reported at steady-state (details in the Supplementary Section 3). Additionally, the isotope fractionation coefficient (ε) was measured at the outlet of the first dilution rate of the chemostat (−5.4‰, only at the lowest dilution rate), and retentostat (−0.45‰).

Calibration Strategy
Our initial intent was to estimate a single set of model parameters for both engineered systems. This was not possible, however, most likely due to differences in bacterial physiology (Ercan et al., 2015;Kundu et al., 2020). In our next attempt, we introduced a switch function (Stolpovsky et al., 2011;Mellage et al., 2015), allowing for environmental-specific transition between the two conditions (chemostat and retentostat) (Supplementary Sections 10.1 and 10.2). This model, despite its high complexity and many degrees of freedom, was still unable to simulate both engineered systems together (Supplementary Section 10.3). Therefore, we exhaustively investigated (using fits for both systems and Frontiers in Environmental Science | www.frontiersin.org September 2021 | Volume 9 | Article 706457 sensitivity analysis) subsets of parameters that could be kept fixed at the chemostat fit while still capturing bacterial behavior in the retentostat in a two-step calibration process, as follows. STEP 1-pre-calibration step: We started by using the five steady-states (one with each dilution rate) measured in the chemostat, and the isotope fractionation of the lowest dilution rate (16 data points). We considered the parameter ranges shown in Table 1. The nominal values were taken from literature ( Table 1). Ranges were selected as to capture parameter variation.
We used the global optimization algorithm Simulated Annealing (simannealbnd) of MATLAB to minimize the weighted sum of squared errors (SSE): where y i obs and y i sim are the mean values per observation type and dilution rate, and the corresponding model output for the i-th data point from n total data points. σ 2 i is the recalculated standard deviation per observation type and dilution rate (details are given in Supplementary Section 3, Supplementary Table 1).
We then calibrated the retentostat system at the steady-state (4 data points) using Simulated Annealing again. An acceptable description could be reached by fixing four parameters and allowing the other four to vary: k AT , K HY M , f cell , and r e (highlighted yellow in Table 1). Details are given in the Section 4.
The model outputs corresponding to the measurements were: Isotope fractionation was determined as: where δ inlet is the isotope ratio of the heavy and the light isotopologues of AT at the inlet, given as −29‰ (Ehrl et al., 2018b, and δ outlet was determined as where R is the reference 13 C/ 12 C isotope ratio of Vienna Pee Dee Belemnite (Brand et al., 2010). The parameter β [Eq.
(2)] can be directly derived from the enzymatic fractionation coefficient of AT (ε −5.4‰) measured for a particular bacterial strain (Ehrl et al., 2018b;Gharasoo et al., 2019): STEP 2-Full calibration: For both systems, a full calibration step, including parameter and output uncertainty were determined with the Markov Chain Monte Carlo (MCMC) algorithm of the DREAM MATLAB toolbox (Vrugt, 2016). We fitted the 8 chemostat system parameters and the 4 differing retentostat system parameters simultaneously (marked in yellow in Table 1) in one optimization run. We chose a flat and uninformative prior distribution for the MCMC. The starting values of the MCMC chains were drawn from a normal distribution of the parameters in log-space (mean value equal to the best fit of the Simulated Annealing (step 1), an arbitrary variance of 1, and zero covariance between the parameters). Minimum and maximum parameter values were taken from Table 1, and the option "reflect" was selected as a method for handling parameter boundaries. TheR-diagnostic (Gelman and Rubin, 1992) lower than 1.2 (Vrugt, 2016) was used as convergence diagnostics. We used a Gaussian likelihood considering heteroscedastic measurement errors as implemented in DREAM:

Exhaustive Soil Extraction
Pesticides (atrazine and hydroxyatrazine) were extracted from soil with an accelerated Solvent Extractor (ASE 300 Dionex, Thermo Scientific) at 80°C and 150 bar, using acetone as the main solvent (parameters in Table 2). To ensure a homogeneous flow through the extraction cells, soil samples were mixed with 80% (mass) clean quartz sand before extraction. To control for potential losses of pesticide during the processing (enrichment and clean-up) of the extracts, 10 ng of Isoproturon-D6 were added to each extract. Subsequently, the extracts were reduced with a rotational evaporator until acetone was evaporated completely. The residual aqueous sample was filtered through 0.25 µm PTFE syringe filters (Agilent, Waldbronn Germany) and 10% (Vol.) of MeOH was added before the measurement at the liquid chromatography-mass spectrometry (HPLC-MS/MS). The target compounds were separated with an Agilent 1290 Infinity HPLC (Agilent, Waldbronn, Germany) using a reversed phase column (Agilent Poroshell 12 EC-C18, 2.7 µm, 2.1 × 100 mm). The quantification of the target compounds was done based on an external calibration using 10 standards with concentrations between 0.02 and 10 μg L −1 . As control for a potential shift during the measurement, every 15 samples, one external standard was monitored, with a concentration of 2.5 μg L −1 (Measurements are shown in Supplementary Section 4, Supplementary Table 2.)

Sorption Test
Six initial concentrations of atrazine (0.06, 0.4, 4, 36, 420 and 2060 μg L −1 ) were prepared from a stock solution of atrazine in MilliQ water (using a pure, analytical standard from Sigma Aldrich). The solutions were spiked with CaCl 2 (0.5 g L −1 ) and NaN 3 (0.25 g L −1 ) to provide a stable ionic strength and minimize bacterial activity. The sorption test was conducted in triplicates in 50 ml glass vials (with Teflon-lined caps), containing 15 g of soil and 30 ml of spiking solution. The vials were kept on a horizontal shaker (150 rpm) for 10 days in the dark and at 20°C. To separate soil solids from water, the vials were kept standing for 3 days until all fine particles were settled. A small test with filtering the aqueous phase had confirmed this approach as valid. Subsequently, the aqueous phase was transferred into clean vials using glass pipettes.
After separating soil solids from water, 20 ng of atrazine-D6 was added as an internal standard to the aqueous phase. Processing of the aqueous samples varied for the different concentrations: Samples with lowest concentrations were enriched via solid phase extraction (Waters OASIS HLB). Samples with expected concentrations between 0.2 and 10 μg L −1 were filtered through 0.25 μm PTFE syringe filters and 2% (Vol.) of acetonitrile was added. For concentrations above 10 μg L −1 , the samples were filtered and then diluted with MilliQ: acetronitrile (98:2) before LC-MS/MS measurements. As quality control, blanks with ultra-pure water, leaching blanks with ultrapure-water and soil, and controls with spiking solution without soil were analyzed in triplicates confirming no relevant loss of atrazine or contamination (Supplementary Section 4, Supplementary Table 3).
We determined the Freundlich sorption parameters (K AT and n AT ) for atrazine at both sites by regressing the sorbed concentration on the solution concentration [Eqs.  Table 4). The sorption coefficient of hydroxyatrazine (K HY ) was calculated by dividing the normalized sorption coefficient of atrazine K * AT (sorption coefficient K AT divided by the water solubility of atrazine S AT ) by the water solubility of hydroxyatrazine (S HY ) at the power of n AT [Eq. (44)] (Carmo et al., 2000;Kleineidam et al., 2002). The sorption exponent for hydroxyatrazine was assumed to be equal to atrazine because the Freundlich exponent is rather soil-than compound-specific:

Soil Predictions
We ran simulations in soils using both sets of calibrated parameters ( Table 4, and Supplementary Section 5, Supplementary  Figure 3). We therefore used the average values for simulating AT and HY fate. We derived the mean water flow in the plough layer (v v ) by dividing the mean daily water flux (0.56 mm d −1 ) by the ploughing depth (300 mm). We fixed the minimum bacterial biomass in soil M according to Klier et al. (2008). The values of the soil parameters are shown in Table 3: To compare with the field monitoring data from the sites Poltringen and Tailfingen, we ran simulations with all four variants of the soil model, assuming an initial application of 1000 μg kg −1 (Vryzas et.al., 2007;Krutz et.al., 2010) and predicting residual concentrations after 30 years.

Global Sensitivity Analysis
We determined the Morris and Sobol indices (Morris, 1991;Campolongo et al., 2007;Saltelli et al., 2008;Pianosi et al., 2015) for the two core model variants (M and T), using the SAFE toolbox of MATLAB (Pianosi et al., 2015;Sigmund et al., 2020). We calculated the mean of the elementary effects (μ*) and the standard deviation of the elementary effects (σ) for the Morris Method, as well as main and total effects for Sobol indices with a total of 15,000 sample inputs in both cases.
We sampled parameters from a uniform distribution taken from the posterior distribution of the fitted parameters against the chemostat and retentostat data combined (Table 4 from the Section 3). We used a Latin hypercube sampling strategy (Marschmann et al., 2019). Additionally for the Morris Method, we calculated the l 2 norm (l 2 μ *2 + σ 2 ) of each parameter (Campolongo et al., 2007;Link et al., 2018;Marschmann et al., 2019) and considered parameters with l 2norm higher than 0.5 as high leverage.
We selected the following outputs: steady state biomass, AT, and HY (extracellular and intracellular), and isotope fractionation ε [Eq. (40)]. We ran the model for 200 days to guarantee steady-state in the simulations.

Local Sensitivity Analysis
We performed a local parametric sensitivity analysis (Ingalls, 2008;Zi, 2011) for the four soil model variants as described above, based on the best fit against the chemostat and retentostat observations. The target outputs were the residual concentration of AT and HY after 30 years. We addressed all kinetics ( Table 1) and soil parameters (Table 3), as well as the initial AT application.

Calibration to Chemostat and Retentostat Data
The two core model variants behave equally in engineered environments, and so we present the results only for Variant T (Results corresponding to Variant M are presented in the Supplementary Section 6, Supplementary Figures 4,5 and Supplementary Table 5). Following a two-step approach, we calibrated the 8 chemostat system parameters and the 4 differing retentostat system parameters simultaneously.

Concentrations
Our simulations were in good agreement with observed data for the chemostat (Figures 2A-C). After the partial re-calibration, we found acceptable agreement for the retentostat system, but with a slightly higher model output uncertainty for the biomass (this was not unexpected, given the relative lack of data for calibration).

Fractionation
Simulations showed agreement with the observed isotope fractionation for both systems, with slightly higher uncertainty for the retentostat ( Figure 2D). Isotope fractionation of AT occurs when enzymatic transformation is the rate-limiting step. In this case, the enzymatic fractionation coefficient of AT (ε) lies close to −5.4‰ (chemostat). At low AT concentrations, the mass transfer across the cell membrane becomes ratelimiting, and no isotope fractionation is observed (ε of just −0.45‰; retentostat) (Ehrl et al., 2018b;Gharasoo et al., 2019;Kundu et al., 2019).

Parameter Estimates and Uncertainty
Kinetic parameters related to AT and HY degradation (chemostat: k AT , K AT M , k HY , K HY M , retentostat: k AT , K HY M ) appear to be well-informed by the data, showing relatively wellconstrained posterior distributions ( Figures 3A,B,C,D,I,J), low standard deviations (Table 4), and considerable impact on model outputs according to the Sobol analyses (especially k AT and K AT M , Supplementary Figure 7). The maintenance parameter m was interestingly well constrained by the chemostat data ( Figure 3E, Table 4); the global sensitivity analysis confirmed this parameter to be low leverage (Supplementary Figures 5,6). The mass-transfer rate parameter r e was not well-constrained for the chemostat data ( Figure 3H), but fitted relatively well to the retentostat data ( Figure 3L), especially with the model variant M (Supplementary Section 6, Supplementary Figure 5). This parameter showed a considerable impact on model outputs (Supplementary Figures 6,7). The yield parameter Y and conversion parameter f cell were highly uncertain and not wellconstrained for either system, probably due to the high correlation with other parameters like kinetic parameters k AT and k HY (Supplementary Tables 6,7 for model variant M).

Comparison of Parameter Estimates Between the Chemostat and Retentostat
Comparing the mean and MAP calibrated parameter values in Table 4, we see that the per-cell AT degradation rate (k AT ) is estimated to be higher for cells living in the low nutrient  Frontiers in Environmental Science | www.frontiersin.org September 2021 | Volume 9 | Article 706457 8 retentostat system. Conjectured physiological adaptations (Kundu et al., 2020) in the retentostat environment may be responsible for the difference in the estimated value of f cell compared to the chemostat, reflecting changes in cell volume. However, this estimate is highly uncertain and highly correlated to other parameters' values. Physiological adaptations might also be responsible for a reduced value of parameter K HY M in the retentostat system, possibly reflecting a change in nutrient demand. The estimate of r e was higher in the chemostat than in the retentostat, indicating a change in membrane properties leading to strong mass transfer limitations across the cell membrane. The estimates of the parameters k AT , K HY M and r e using model variant M show the same tendencies, but exhibit stronger changes (increase/decrease) from chemostat to retentostat (Supplementary Section 6, Supplementary Table 5). The main difference is in parameter f cell , which shows a clear reduction in the retentostat, strongly supporting the findings of Kundu et al. (2020).

Predictions of Atrazine and Hydroxyatrazine Fate in Soils and Comparison Against Field Data
We simulated the fate of AT and HY in soil for 30 years, assuming a single initial AT input of 1,000 μg kg −1 (Vryzas et al., 2007;Krutz et al., 2010). For this, we used the full posterior parameter estimates from the chemostat and retentostat systems for four model variants (Figure 4). All model variants predicted very low residual AT concentrations, considerably underestimating the observed concentrations of 0.3 and 0.6 μg kg −1 in the top soil of both field sites (Poltringen and Taiflingen respectively) ( Figures 4A,B).
In contrast, predictions of residual HY mainly overestimated the observed HY concentration at both study sites (around 2 μg kg −1 in both sites) ( Figures 4C,D). Predictions using retentostatfitted parameters in combination with thermodynamically constrained growth and leaching ( Figure 4D) predicted longterm persistence of HY, with mean values around 36 μg kg −1 . However, model variants with Monod kinetics (M-NL(R) and M-L(R)) performed better and predicted residual HY concentrations much closer to the measurements (9 and 20 μg kg −1 , Figure 4D).
As expected, simulations of this simple model over 30 years are highly uncertain. Based on our local sensitivity analysis (Supplementary Section 9,Supplementary Tables 8,9), the sorption exponent of both chemicals (n AT and n HY ) showed the highest impact on the residual concentrations of AT and HY after 30 years, revealing a strong dependency on sorption characteristics of the soils. Surprisingly, the initial application of AT only impacted the residual concentration of AT in model variants incorporating thermodynamic growth constraints. Water flow (v v ) and minimum bacterial biomass (M) had low impact on the residual concentrations, despite their role to improve model predictions (our best model predictions  include leaching; recall that parameter M accounts for alternative carbon sources for soil bacterial biomass). As to be expected, the kinetic parameters, in contrast to the sorption parameters, had a negligible impact on the target outputs.

Bacterial Adaption to Low Nutrient Availability Affects Model Parameterization
Due to the apparent similarities between the chemostat and retentostat systems, our initial intent was to achieve a joint fit for both systems. In particular, by including a flexible formulation of the mass-transfer rate, as well as a thermodynamically constrained growth rate instead of a Monod formulation, we aimed to represent systems with or without mass-transfer limitations across the cell membrane by one model. However, we found that goal unattainable. Recent publications (Ercan et al., 2015;Kundu et al., 2020) show evidence of a phenotypic differentiation of a single population into separate growing and non-growing (i.e. energy used only for maintenance) bacterial subpopulations (Kundu et al., 2020). Thus, we focused on the key parameters that have to be re-calibrated between the two systems using two model variants that exhibit equivalent behavior over the range of inputs in the engineered systems (Table 4).
After fitting model parameters to the chemostat data, we systematically tested which parameters had to be re-calibrated to capture the retentostat behavior. We were guided by sensitivity analyses, as well as our understanding of the role of the parameters in our model. We fixed the maximum degradation rate of HY (k HY ), the growth yield (Y), the half-saturation concentration for AT degradation (K AT M ) and the maintenance parameter (m) because of their low impact on model output (Supplementary Material, Supplementary Figures 5,6). Similar sensitivities were previously reported in the literature . Summing up, the parameters that had to be recalibrated to capture the retentostat behavior are: k AT , K HY M , f cell , and r e . We justify the requirement of these needed adjustments in the following.
The parameters k AT (maximum degradation rate of AT), and K HY M (reference/half-saturation concentration) represent physiological features that can be expected to change under starvation conditions (Lever et al., 2015;Kundu et al., 2020). Relative to the chemostat conditions, in the low-HY retentostat environment, we estimate a higher values of k AT and lower values of K HY M ( Table 4), indicating faster AT transformation to HY, and physiological adaptation of microorganisms to use of HY, respectively. While the fitted value of k AT was about twice the value of k HY (maximum degradation rate coefficient of HY) in the retentostat, both parameters (k AT and k HY ) were similar in the chemostats ( Table 4). This difference in the parameterization of both systems shows that the physiological adaption of microorganisms to low concentrations affects the regulation of the AT degradation reaction network such that HY transformation becomes rate-limiting for microbial growth.
We found that re-calibration of the parameter f cell is an efficient way to capture specific bacterial differentiation for low nutrient systems . The parameter f cell is a scaling factor used to convert cells to C (Vrede et al., 2002) and might suggest morphological changes (shape and volume) observed in Arthrobacter aurescens to cope with stressful starvation conditions (Mongodin et al., 2006;Lever et al., 2015). Due to the high uncertainty in parameter estimation, more experiments are needed to identify the underlying mechanism.
Changes in the value of r e (mass-transfer rate coefficient) between chemostat and retentostat system could reflect morphological/physiological changes in the cell membrane ( Table 4). The relatively lower value of r e in the retentostat suggests a strong mass-transfer limitation across the cell membrane in that case.

Pesticide Persistence in Soil
The main objective of our work was to accurately represent Atrazine (AT) degradation in soils, and especially to capture the long-term persistence of AT and its main metabolite Hydroxyatrazine (HY).
Despite the related uncertainty for long-term predictions, persistence of HY even after 30 years was consistently predicted by model variant M-L calibrated with retentostat data ( Figure 4D). In general, retentostat concentrations are closer to the soil environment, so that more accurate predictions are to be expected (biomass retention, low nutrient levels). Additionally, incorporation of leaching gave a better representation of the pesticide losses over time. Simulation with a simple model incorporating only leaching over the 30 years leads to a residual concentration of AT of about 2 μg kg -1 . This value is close to the measured residual concentrations indicating that only low AT degradation might have occurred at the field sites ( Figure 4). Standard Monod model variants predicted HY concentrations after 30 years better than thermodynamic models. Therefore, energetic constraints of microbial growth likely do not limit HY degradation in soils.
In contrast, all model configurations predicted a nearly complete consumption of AT after 30 years, a behavior not observed in field surveys (Jablonowski and Schäffer, 2011;Vonberg et al., 2014), including the field measurements of this study ( Figure 4A).
A range of biological and physical processes in soil have been hypothesized as potential mechanisms of pesticide persistence in natural systems. These include physicochemical control mechanisms limiting bioavailability, such as chemisorption onto humic substances (Mudhoo and Garg, 2011), physical stabilization in soil micro aggregates (Hochman et al., 2021), or the spatial encounter of substrates and degraders (Shi et al., 2021). Including these additional mechanisms by applying better sorption and stabilization model formulations (van Genuchten and Wagenet, 1989;Streck et al., 1995;Altfelder and Streck, 2006;Villaverde et al., 2009;Kästner et al., 2014;Yu et al., 2020) and spatially resolved modeling approaches (Babey et al., 2017;König et al., 2020;Pagel et al., 2020;Pot et al., 2021) might further improve predicting the persistence of AT and other pesticides in soil. Our study investigated to what extent mass-transfer limitations and bioenergetic constraints can explain the longterm fate of atrazine and its major metabolite hydroxyatrazine in soils. We found evidence against the hypothesis that passive diffusion across the cell membrane of bacterial degraders limits atrazine degradation in the long term. Atrazine is not degraded to HY for the energy gain by microorganisms and our results suggest that sorption-limited bioavailabilty and not energetic growth constraints control the persistence of hydroxyatrazine. Hence, standard Monod kinetics for bacterial growth can predict the long-term fate of organic chemicals if soil microorganisms directly utilize them as an energy source. Further research should prioritize the analysis of energetic costs of biogeochemical transformations without a direct microbial energy gain (atrazine dechlorination).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
LR, HP, TS, and BI contributed the conception of the manuscript and the design of the study. LR developed the model, performed the simulations, and analyzed the data. KK provided the chemostat/retentostat data. JM performed soil and seepage water measurements (sampling, processing, and analysis). LR wrote the first draft of the manuscript. All authors contributed to the manuscript revision, read and approved the submitted version.