Calculation and Interpretation of Substrate Assimilation Rates in Microbial Cells Based on Isotopic Composition Data Obtained by nanoSIMS

Stable isotope probing (SIP) combined with nano-scale secondary ion mass spectrometry (nanoSIMS) is a powerful approach to quantify assimilation rates of elements such as C and N into individual microbial cells. Here, we use mathematical modeling to investigate how the derived rate estimates depend on the model used to describe substrate assimilation by a cell during a SIP incubation. We show that the most commonly used model, which is based on the simplifying assumptions of linearly increasing biomass of individual cells over time and no cell division, can yield underestimated assimilation rates when compared to rates derived from a model that accounts for cell division. This difference occurs because the isotopic labeling of a dividing cell increases more rapidly over time compared to a non-dividing cell and becomes more pronounced as the labeling increases above a threshold value that depends on the cell cycle stage of the measured cell. Based on the modeling results, we present formulae for estimating assimilation rates in cells and discuss their underlying assumptions, conditions of applicability, and implications for the interpretation of intercellular variability in assimilation rates derived from nanoSIMS data, including the impacts of storage inclusion metabolism. We offer the formulae as a Matlab script to facilitate rapid data evaluation by nanoSIMS users.


INTRODUCTION
Stable isotope probing (SIP) is an experimental approach used to trace the fates of substrates within a microbial community. In this approach, samples are incubated with a target substrate that is enriched in a stable isotope, such as 13 C, and then cells or cellular components are analyzed for enrichment in the target isotope (Boschker and Middelburg, 2002;Dumont and Murrell, 2005;Neufeld et al., 2007;Hungate et al., 2015). SIP can be used to identify microbes in environmental samples that are actively using a particular growth substrate by selectively recovering and analyzing isotopic signatures in biomarkers such as DNA, rRNA, lipids, fatty acids or proteins (Boschker and Middelburg, 2002;Neufeld et al., 2007;Jehmlich et al., 2016). SIP can also be used to quantify rates of substrate assimilation and compartmentalization by specific community members with the goal of estimating fluxes at the ecosystem scale (Dumont and Murrell, 2005;Hungate et al., 2015). In the latter application, use of nano-scale secondary ion mass spectrometry (nanoSIMS; see Hoppe et al., 2013or Nuñez et al., 2017 for a review of the method) to measure substrate assimilation into individual cells has become popular in the last decade (Musat et al., 2012;Pett-Ridge and Weber, 2012;Mayali, 2020;Ploug, 2021). Coupling SIP with nanoSIMS enables measurement of activities of microbes from the environment, including those that cannot yet be cultured (Krupke et al., 2015;Harding et al., 2018;Dekas et al., 2019;Geerlings et al., 2020;Mills et al., 2020), are rare (Musat et al., 2008;Zimmermann et al., 2015), or are tightly associated with their geophysical and/or biological environments (Foster et al., 2011(Foster et al., , 2013Milucka et al., 2012;Bonnet et al., 2016;Berthelot et al., 2019;Loussert-Fonta et al., 2020), thus offering opportunities to advance understanding of how single-celled life works, adapts, and impacts geochemical element cycling.
A typical output of a nanoSIMS measurement is the isotopic composition of a cell, determined as the cell-specific isotope ratio, R, or atom fraction, x [e.g., for carbon isotopes, R = 13 C/ 12 C and x = 13 C/( 12 C+ 13 C)]. One way of using this data is to convert R or x into quantities that can directly be interpreted in terms of the cell's mass balance, including Fx net , K A or X net (Popa et al., 2007;Finzi-Hart et al., 2009;Stryhanyuk et al., 2018;Dekas et al., 2019). Specifically, the parameter Fx net , introduced by Popa et al. (2007) and equal to the parameter K A introduced by Stryhanyuk et al. (2018), reflects the net amount of element (e.g., carbon) assimilated by the cell during the SIP experiment (E a ) relative to the initial content of the element in the cell prior to the SIP experiment (E i ), i.e., Fx net = K A = E a /E i . Similarly, the parameter X net , used for instance by Dekas et al. (2019), reflects the net amount of element assimilated by the cell relative to the final element content in the cell (E f ), i.e., X net = E a /E f ( Table 1).
The quantities described above do not account for incubation duration so can only be compared across SIP experiments that are incubated for the same duration. Calculating the rate of substrate assimilation obviates this requirement and has become a common output from nanoSIMS analyses of microbial cells (Ploug, 2021, and references therein). The two types of substrate assimilation rates used are the substrate-specific and cell-specific assimilation rates. While the substrate-specific rate [e.g., the rate of carbon assimilation normalized to the carbon content of the cell; in mol C (mol C) −1 h −1 ] provides useful information on cellular turnover of the substrate and can be converted to a cell's doubling time (e.g., Musat et al., 2008;Martínez-Pérez et al., 2016;Eichner et al., 2017), the cell-specific rate (e.g., the carbon assimilation rate per cell; in mol C cell −1 h −1 ) can be used to estimate the impact of a cell's population on element fluxes over larger spatial and temporal scales (e.g., Foster et al., 2011;Krupke et al., 2015;Bonnet et al., 2016;Klawonn et al., 2016;Arandia-Gorostidi et al., 2017;Mills et al., 2020).
Calculating substrate assimilation rates by individual cells from nanoSIMS data requires consideration of cell division and the mathematical description of cell growth. Substrate assimilation during a SIP incubation, as indicated by isotopic enrichment of cell biomass at the end of the incubation (i.e., a positive value of the cell-specific Fx net or X net ), implies that there is a non-zero probability that the cell divided during the SIP incubation. A simple mass balance calculation reveals that, under balanced growth, this probability reaches 1 when the cell-specific quantities Fx net and X net exceed 1 and 0.5, respectively. Some recent studies have considered cell division occurring during SIP incubations by applying a correction factor (i.e., 2) in growth rate calculations (Schoffelen et al., 2018(Schoffelen et al., , 2019Olofsson et al., 2019a,b). The mathematical basis for such correction factors has not been, however, explicitly formulated, and the impact of cell division on calculated rates of substrate assimilation in cells has not been investigated in a systematic way.
Calculated substrate assimilation rates will vary depending on the model used to describe substrate assimilation over the cell cycle (i.e., cell growth) and thus over the incubation time. A linear model assumes that substrate assimilation proceeds at a constant cell-specific rate (i.e., zero-order kinetics) leading to a linear increase in the biomass of individual cells over time. In contrast, an exponential model assumes that substrate assimilation proceeds at a rate that is linearly proportional to the instantaneous cell biomass (i.e., first-order kinetics) leading to an exponential increase in the biomass of individual cells over time (Collins and Richmond, 1962;Koch and Schaechter, 1962;Koch, 1966). While the linear model of substrate assimilation is most commonly used to calculate assimilation rates in cells from nanoSIMS data (e.g., Foster et al., 2013;Klawonn et al., 2016;Stryhanyuk et al., 2018;Mills et al., 2020), the exponential model is used rather rarely (e.g., Martínez-Pérez et al., 2016;Berthelot et al., 2019;Geerlings et al., 2020Geerlings et al., , 2021Polerecky et al., 2021). The impact of the substrate assimilation model (linear vs. exponential) on the calculated assimilation rate has previously not been studied in detail.
Assimilation rates determined at single-cell resolution often reveal intercellular heterogeneity in populations. Previous nanoSIMS studies typically attributed such heterogeneity to factors that affect the assimilation of an (isotopically labeled) element during the incubation (causing variation in E a ), i.e., differences in the intrinsic metabolic activity of a cell, which may vary with cell cycle stage or cell age, or stochastic gene expression (e.g., reviewed by Ackermann, 2015). Other sources of heterogeneity have received less attention, such as assimilation of unlabeled sources of the target element, which will affect the degree of isotopic labeling of the assimilated material, or differences in the cell's initial elemental content (variation in E i ). For instance, selective assimilation of C and N into storage inclusions or mobilization of C or N from existing cell material during a SIP experiment can critically affect both the (apparent) amount of an element assimilated during an incubation (E a ) and the initial amount of an element present in the cell (E i ) (Polerecky et al., 2021). Thus, metabolism of storage compounds can partly

Symbol
Unit Description C, C i mol C cell −1 Carbon content of a cell; 'i' refers to the initiation of the SIP incubation, i.e., time-point t = 0.
C max mol C cell −1 Carbon content of a cell just before binary cell division; C max /2 corresponds to the carbon content just after binary cell division.
s -Cell cycle stage, calculated from the C content as s = C/(C max /2) -1.
C mol C cell −1 Carbon content of a cell averaged over the cell cycle. It is equal to the average in a population of cells with perfectly unsynchronized cell cycles. It is related to C max according to C = C max /(2 · ln (2)) and C = C max · ln (2) for cells assimilating C with zero-order and first-order kinetics, respectively (Koch, 1966).
x E S -Source-normalized excess 13 C atom fraction of a cell (Eq. 10).
r mol C cell −1 h −1 Cell-specific rate of carbon assimilation by a cell.
Average cell-specific C assimilation rate in a population.
k h −1 Carbon-specific rate of carbon assimilation by a cell. t h Doubling time; time needed for a cell to double its carbon content. Calculated as τ = C max /(2 · r) and τ = (ln2)/k for a cell assimilating C with zero-order and first-order kinetics, respectively.
ρ mol C µm −3 Carbon density of a cell.
x S,tar -13 C atom fraction of the target source provided externally during the SIP incubation.
x S,alt -13 C atom fraction of an alternative carbon source; can be external (e.g., present in the environment of the cell) or internal (present in the form of intra-cellular C storage inclusions).
x S,eff -Effective 13 C atom fraction of the carbon source (Eq. 6).
f tar -Fraction of C assimilated by the cell from the target source; the remaining fraction of C is assimilated from the alternative C source (f alt = 1f tar ).
θ -Probability density function (PDF) describing the distribution of C among cells in a population. Given by Eq. 16 and 18 for a population with perfectly unsynchronized and partially synchronized cell cycles, respectively. Examples shown in Figure 4B, Supplementary Figures 3 and 4).
ζ -Probability density function (PDF) describing the distribution of x E S in a population of cells assimilating C with zero-order kinetics. Examples shown in Figure 4C. explain variability in the enrichment of cellular biomass (e.g., Fx net or X net ). In this work we use mathematical modeling to identify key processes and parameters that play a role when calculating and interpreting cellular substrate assimilation rates from the isotopic composition information determined by SIP-nanoSIMS. First, we investigate how the isotopic composition of a cell will vary in time depending on the model of substrate assimilation (zeroorder vs. first-order kinetics), the substrate assimilation rate, cell biomass, and the isotopic composition of the substrate. We specifically analyze how cell division during the incubation impacts the predicted isotopic composition of a cell and how the isotopic compositions will vary among cells within a population depending on the degree of cell cycle synchronicity. Based on these theoretical analyses, we provide nanoSIMS users with step-by-step guidelines for calculating substrate assimilation rates including formulas and their underlying assumptions. We also identify key biological factors that can introduce variability into the calculated rates, specifically highlighting the previously unrecognized role of storage inclusions and their impact on interpretations of intercellular variability in calculated assimilation rates.

MATERIALS AND METHODS
Our analysis is based on mass balances, which we formulate as differential equations to facilitate identification of parameters impacting the predicted isotopic composition of cells during a SIP experiment. We focus our analysis on assimilation of carbon (C) from a 13 C-labeled carbon source. Analogous analyses can be conducted for other elements and their stable isotopes commonly used in SIP experiments (e.g., 15 N, 18 O). Throughout this study, we generally adhere to the accepted notation guidelines for reporting isotope enrichment data (Coplen, 2011) with a few exceptions (see Table 1 for a list of symbols and their definition used in this study).

C Labeling of an Individual Cell
The model of C assimilation by an individual cell is mathematically formulated by a differential equation where C denotes the cell-specific C content (in mol C cell −1 ) and r denotes the cell-specific C assimilation rate (in mol C cell −1 h −1 ). Both quantities are a function of time, as indicated by the argument (t). In this study, we neglect the potential effect of kinetic isotopic fractionation as it is typically minor compared to the changes in the isotopic composition of the cell due to labeled substrate assimilation during the SIP incubation. If needed, the effect can easily be included in the analysis described below. Additionally, we do not consider processes associated with cellular biomass turnover, e.g., exudation of freshly fixed C or respiration of storage compounds. Such processes could be implemented by assuming a second, "loss" term in Eq. 1. Their mathematical treatment is, however, beyond the scope of this study. Given the rate of increase in the total C (i.e., 12 C + 13 C) content of the cell, r, the 13 C content of the cell will increase at a rate r multiplied by the probability that the assimilated C atom is 13 C. This probability is equal to the 13 C atom fraction of the carbon source, denoted as x S . Consequently, the 13 C content of the cell is described by a differential equation To make our analysis more general, we assume that the assimilated C may originate from two sources: (1) the target C source, which is isotopically labeled by the experimenter by the addition of a naturally less abundant stable isotope (e.g., 13 C), and (2) an alternative C source, which only contains the natural abundance of 13 C and may be present and assimilated by the cell without experimenter's knowledge. For example, 13 Clabeled dissolved inorganic carbon (DIC) is a common target C source and unlabeled dissolved organic carbon (DOC) could be an alternative C source. Assuming that C assimilations from the target and alternative sources are independent of each other, the total C and 13 C contents of the cell will vary with time according to differential equations where subscripts "tar" and "alt" refer to the target and alternative C source, respectively. Although in some SIP experiments the isotopic composition of the C sources may vary in time (e.g., Schreiber et al., 2016;Geerlings et al., 2020), we consider hereafter the more common situation that x S,tar and x S,alt are time-independent. We make this assumption because our aim is to derive simple formulas to calculate the substrate assimilation rate based on the isotopic composition of the cell. A more advanced analysis describing SIP experiments with time-dependent substrate labeling is provided in Supplementary Methods, Section "Modeling Cellular Assimilation of Substrates With Time-Dependent Isotopic Composition, " and in Supplementary  Figures 1,2. Considering this assumption, we denote the total cell-specific C assimilation rate as r = r tar + r alt , and the fraction of C assimilated by the cell from the target source as f tar = r tar /(r tar + r alt ), which makes Eq. 3 identical to Eq. 1 and yields the following differential equation for the cellular 13 C content: where is the effective 13 C atom fraction of the carbon source.
To understand the consequences of C assimilation from a source characterized by x S,eff for the isotopic composition of the cell requires the solutions of differential equations 1 and 5. The result depends on how r may or may not vary across the cell cycle. We consider two idealized scenarios: (1) r is constant across the cell cycle (i.e., the cell assimilates C according to zero-order kinetics), (2) r is linearly proportional to the instantaneous cellspecific C content, r = k · C (i.e., the cell assimilates C according to first-order kinetics), where k is constant across the cell cycle. The model of zero-order kinetics can be conceptualized as a cell whose C content increases due to the constant activity of the cell components (e.g., ribosomes, proteins) that are present at the beginning of the cell cycle and do not change in concentration over its course ( Figure 1A). The model of first-order kinetics is conceptualized as described by Koch and Schaechter (1962) where the cell components are synthesized continuously and operate at constant rates throughout the cell cycle such that the ratio between these operational units and the mass of the total cell remains constant across the cell cycle ( Figure 1B). Additionally, we assume that the C density of a cell, denoted hereafter as ρ (in mol C µm −3 ), is constant across the cell cycle. Variation in ρ will occur if C is assimilated into storage inclusions (e.g., polysaccharide or lipid inclusions) or if C storage inclusions FIGURE 1 | Conceptual diagram of C assimilation during cell growth. Cell is represented by a rectangle with round edges. Size of the rectangle represents C content of the cell, thickness of the arrow entering the cell represents the cell-specific rate of C assimilation, r. (A) Zero-order kinetics of C assimilation assumes that r is independent of the cellular C content. (B) First-order kinetics of C assimilation assumes that r increases proportionally to the instantaneous cellular C content.
Frontiers in Microbiology | www.frontiersin.org are catabolized for respiration or the synthesis of other cell components during certain periods of the cell cycle. Effects of these variations are qualitatively examined in the Discussion.

Zero-Order Kinetics of C Assimilation
Assuming constant r, solutions to the differential equations 1 and 5 are readily found as where C i and x i denote, respectively, the C content and 13 C atom fraction of the cell at the time-point when the SIP experiment with the labeled substrate was initiated (i.e., at t = 0). Thus, zeroorder kinetics of C assimilation corresponds to a linear increase in cell biomass.
Equations 7-8 imply that the 13 C atom fraction of the cell, defined as x = 13 C/( 12 C+ 13 C), increases in time according to To simplify this function, we define the source-normalized excess 13 C atom fraction as This is a convenient quantity for describing the isotopic composition of the cell because it increases from 0 to 1 as the cell increases its 13 C atom fraction from x i towards x S,eff (see Supplementary Methods, Section "Relationships between x E S , X net , Fx net and K A , " for a relationship between x E S and the quantities Fx net , K A and X net introduced previously by other authors). Upon rearrangement, this definition combined with Eq. 9 yields The next step is to account for cell division. Specifically, we assume binary division, in which the cell divides into two identical cells when its C content reaches some critical value C max (Collins and Richmond, 1962;Koch, 1966;Westfall and Levin, 2017). Hence, an important assumption underlying the models in this study is that cells are not starving, a condition that is commonly associated with accumulation of storage products with no regular cell division. Thus, equations 7-11 correctly describe the C content and 13 C atom fraction of a cell if, and only if, C(t) is in the interval between C max /2 and C max . When C(t) in equations 7-11 reaches C max , it must be reset to C max /2. As a result, 13 C labeling of the cell, described by the parameter x E S , will vary in time following a "zig-zag" pattern that depends on C i , C max , and r (see Results, Figure 2A). Hereafter, we denote the function describing this pattern as Z(t, C i , C max , r). The consequences of the zig-zag pattern on the resultant substrate assimilation rate are discussed in Results. Note that although zero-order kinetics of C assimilation implies a linear increase in the biomass of an individual cell (Eq. 7), binary cell division of such cells yields a population growing exponentially in cell number and total cell biomass over time.

First-Order Kinetics of C Assimilation
Under first-order kinetics, r is linearly proportional to the instantaneous cell-specific C content, i.e., r (t) = k · C (t), where k is a constant. The parameter k is commonly referred to as the growth rate constant (in h −1 ) and represents the carbon-specific rate of C assimilation [in mol C (mol C) −1 h −1 ]. Substituting this rate expression into Eq. 1, we obtain for the cell-specific C content where C i is the initial C content of the cell. Thus, firstorder kinetics of C assimilation corresponds to an exponentially growing cell. Note that in some studies (e.g., Martínez-Pérez et al., 2016) exponential cell growth is assumed to be described by a function 2 Growthrate·t rather than e k·t . These descriptions are equivalent when "Growthrate" is calculated as k/ln (2). The corresponding doubling time, τ, is calculated by considering that C(τ) = 2·C i , which yields a formula τ = ln(2)/k = 1/Growthrate. To evaluate the dynamic of the 13 C atom fraction (x) for an exponentially growing cell, we start from the definition of x and apply the quotient rule to evaluate dx/dt. Then, we use equations 1, 5 and 12 yielding the following differential equation for x: The solution to this differential equation is Using the definition in Eq. 10, this function is simplified to when the isotopic composition of the cell is described by x E S instead of x. Thus, 13 C labeling of a cell assimilating C according to first-order kinetics is described by one minus an exponential function of time.
To account for cell division, we again assume that the cell divides into two identical cells when its C content reaches some critical value C max . Thus, when C(t) described by Eq. 12 reaches C max , it must be reset to C max /2. In contrast to cell-specific C assimilation modeled by zero-order kinetics, this reset in C(t) due to cell division does not affect the time-dependence of the isotopic composition of the cell assimilating C according to first-order kinetics. This is because the assimilation rate decreases by half as the cell divides, hence the ratio of 13 C uptake to 12 C present in the cell remains the same. This insight also emerges directly from the differential equation 13, which shows that x, and hence also x E S , is independent of C i (see Eq. 14-15).

C Labeling in a Population of Cells
Here, we expand our analysis from the single cell view to a population of cells with varying degrees of cell cycle synchronization. Our aim is to reveal how the predicted isotopic FIGURE 2 | Simulations of the time-dependence of the isotopic composition and C content in cells assimilating C according to different models. Shown are the source-normalized excess 13 C atom fraction (panels A,D), total C content (panels B,E), and 13 C content (panels C,F) for two cells with different initial C contents (C i = C max /2 for cell 1, C i = 1.5·C max /2 for cell 2). Note that cells divide when their total C content reaches C max . For comparison, results calculated based on a model that ignores cell division are shown by dashed lines in panel (A). The number of primes behind the digit 1 and 2 indicates the generation of the daughter cells corresponding to the original mother cells. The cells assimilating C with zero-order kinetics had equal cell-specific C assimilation rates, r, while the cells assimilating C with first-order kinetics had equal carbon-specific C assimilation rates, k. These rates were related as k = r/ C , where C was the same for both cell types. In all panels the time was normalized by the doubling time τ, calculated as τ = C max / (2 · r) and τ = ln (2) /k for the cells assimilating C with zero-order and first-order kinetics, respectively. Blue arrows with numbers at t = n · τ indicate x E S = 1 − (1/2) n , where n = 1, 2, etc. Simulations were performed for x i = 0.011 and for hypothetical values of x S,eff = 0.2, C = 10 fmol C cell −1 and r = 1.35 fmol C cell −1 h −1 . These values correspond to k = 0.135 h −1 , τ ≈ 5.1 h, and C max ≈ 13.9 and 14.4 fmol C cell −1 for the cells assimilating C with zero-order and first-order kinetics, respectively. enrichment varies within a population of cells assimilating C at equal rates and how this variability depends on the model used to describe C assimilation by the cell and the synchronization of the cell cycles among cells. Note that regardless of whether C assimilation is modeled by zero-or first-order kinetics, binary division at the end of each cell cycle causes the population to grow exponentially. Moreover, the population doubling time, τ, is related to the C assimilation rate according to τ = C max /(2 · r) and τ = ln 2 /k for the models assuming zero-order and first-order kinetics of C assimilation, respectively.

Zero-Order Kinetics of C Assimilation in a Population
First, we consider zero-order kinetics of C assimilation by individual cells. If cell cycles in a population are not synchronized, C at a given point in time varies among cells between C max /2 and C max (Collins and Richmond, 1962;Koch, 1966). Thus, because of the dependence of the Z function on C i (see above, Section "Zero-order Kinetics of C Assimilation"), individual cells will have different x E S depending on their C i even if they assimilate C at the same rate r. We denote the probability density function (PDF) describing the distribution of C i and x E S among cells by θ and ζ , respectively, and first consider perfectly unsynchronized and then partially synchronized cells. In both cases we assume that the critical C content where each cell divides, C max , is equal for all cells.
As shown by Koch (1966), the PDF for perfectly unsynchronized cells is for m in the interval C max /2 ≤ m ≤ C max , and it is equal to zero for m outside this interval. For these cells, the C content of the average cell in the population, which is equivalent to the C content averaged over the cell cycle, is related to the critical C content C max by the following formula (Koch, 1966): Furthermore, both θ (m) and C are time independent as the population grows. For partially synchronized cells, θ can take various forms. Here, we assume θ to be based on a Gaussian function. However, because θ is non-zero within the interval between C max /2 and C max , the tails of the Gaussian function reaching outside this interval must be constrained within this interval with a factor that accounts for binary cell division (Supplementary Figure 3). This consideration yields for m in the interval C max /2 ≤ m ≤ C max , and zero for m outside this interval. In this function, C 0 and C describes the center and width of the Gaussian function, respectively (C 0 must lie between C max /2 and C max ), and A is the normalization constant such that C max C max /2 θ (m) dm = 1. We define the degree of cell cycle synchronicity in the population by the ratio γ = C max /(8 · C) (Supplementary Figure 3). Thus, a greater degree of synchronicity corresponds to a smaller C, i.e., a narrower distribution of C among cells (Supplementary Figures 4A,B), and vice versa. Specifically, at the limits of C → 0 or C → ∞, the function θ(m) in Eq. 18 describes populations of cells with perfectly synchronized (γ → ∞) or perfectly unsynchronized (γ → 0) cell cycles, respectively. In the latter case, θ(m) in Eq. 18 and Eq. 16 become equivalent as required (Supplementary Figure 4C). Note that if γ is large (i.e., C C max ) and each cell in the population grows at the same cell-specific rate r, Eq. 18 implies that 95% of cells divide within a time interval t 95 = 4 · C/r (Supplementary Figure 3). Thus, for a population with cell cycles that are highly synchronized, the parameter γ is equal to the ratio between the population doubling time, τ = C max /(2 · r) ( Table 1), and the interval during which 95% of cells undergo division (i.e., γ = τ/ t 95 ; Supplementary Figure 3).
We performed Monte-Carlo simulations in Matlab to evaluate ζ from θ. Specifically, a cell j with a random initial C content C ij is selected based on the PDF θ. Then, for a given value of r and t, the isotopic composition of daughter cells originating from cell j is calculated as x E Sj = Z t, C ij , C max , r . Additionally, the number of daughter cells originating from cell j is calculated from the number of cell divisions as N j = 2 n j , where n j denotes the largest integer that is smaller than n j = C ij + r · t /(C max /2) − 1. Based on the values of x E Sj and N j obtained for many random choices of C ij , a histogram approximating the PDF ζ is reconstructed (see Results, Figure 4C). Finally, the average x E S for all cells in the population is calculated as where the index k refers to the k th bin in the histogram approximating ζ . Note that both the function ζ and the average value x E S depend on time, t, although this dependence was omitted in this expression to simplify notation.

First-Order Kinetics of C Assimilation in a Population
Equations 14-15 show that isotopic composition of a cell assimilating C according to first-order kinetics is independent of C i . Thus, in contrast to a cell assimilating C according to zero-order kinetics, the average x E S for the population of cells assimilating C by first-order kinetics is described by the same equation as for an individual cell (Eq. 15), and the variance among cells is zero.

Modeling 13 C Labeling in Cells
Equations 9-11 imply that 13 C labeling of a cell assimilating C according to zero-order kinetics depends on the cell-specific rate of C assimilation, r, and the initial C content of the cell, C i . We illustrate the consequences of these dependencies by considering two cells with identical r but different C i . We assume C i1 = C max /2 for cell 1, which corresponds to a cell immediately after division (initial cell cycle stage s i1 = 0), and C i2 = 1.5 · C max /2 for cell 2, which corresponds to a cell in the middle of its cell cycle (s i2 = 0.5; see Table 1 for the definition of s). Both cells then assimilate C at the same and constant rate r over three cell cycles (Figures 2A-C). The lower C i1 compared to C i2 causes x E S1 to initially increase faster than x E S2 (Figure 2A, t/τ = 0.5) and the C content of cell 2 to reach C max earlier than that of cell 1. At C max , the rate of increase in x E S in the two daughter cells of cell 2 abruptly increases 2-fold due to the abrupt decrease in their C content from C max to C max /2 and an unchanged r (Figures 2A,B, t/τ = 0.5). A similar "kink" in the evolution of x E S is observed after cell 1 divides, but it occurs later because of the lower initial C content of cell 1 (Figure 2A, t/τ = 1).
The x E S values for the two cells described in Figure 2 become equal at regular time intervals. These intervals are separated by τ = C max /(2 · r), i.e., the time needed for each cell to double its C content (Figure 2A). This result is expected because the total C content of the daughter cells at time-points t = n · τ (n = 1, 2, . . .) will be the same as the total C content of the original mother cell at t = 0 (Figures 2B,C). Hence, the zero-order kinetics model of C assimilation yields x E S = 1 − (1/2) n at time-points t = n · τ irrespective of the value of C i (Figure 2A).
The time-dependence of x E S in a cell assimilating C according to first-order kinetics is described by a function that is independent of C i (Eq. 15). Setting the initial C contents of two cells to C max /2 and 1.5 · C max /2 leads to different timedependences of the total C and 13 C contents for each cell (Figures 2E,F). However, the increase in x E S over time is the same for both cells (Figure 2D). Thus, cell division and the cell's initial C content do not influence how the cell's isotopic composition varies in time if C assimilation can be adequately captured by the first-order kinetics model. While x E S of an individual cell modeled by zero-order kinetics is described by a function with abruptly changing time-derivative (Figures 2A, 3A), x E S for the average cell representing its population is described by a smooth function. If cell cycles in the population are perfectly unsynchronized, the function is identical to the function describing x E S of a cell modeled by first-order kinetics with k = r/ C , where C is the average C content of the cells across the cell cycle (compare solid and dashed lines in Figure 3B). Thus, for an average cell in a population with perfectly unsynchronized cell cycles, both the C content and 13 C labeling behave in the same way regardless of whether their C assimilation is described by zero-order or firstorder kinetics. This important result underpins the procedure for calculating rates of C assimilation in cells based on nanoSIMS data (see next section). If cell cycles in a population are partially synchronized (as is typical, e.g., for photoautotrophs that tune cell division according to day-night cycles), the average x E S predicted by the zero-order kinetics model follows a zig-zag pattern like that describing x E S for a single cell ( Figure 4A). This behavior is a result of the variation of the probability density function ζ becoming increasingly pronounced with the increasing degree of cell cycle synchronicity ( Figure 4C). However, the interval of x E S where ζ is non-zero is independent of the degree of cell cycle synchronicity ( Figure 4C). Thus, the value of x E S for any cell from a population with partially synchronized cell cycles, including the cell representing the average of the population, always lies within the interval defined by the minimum and maximum value of x E S determined for cells with perfectly unsynchronized cell cycles (Figures 4A, 4C). Note that the width of this interval varies in time following a wax-and-wane pattern oscillating between zero and a maximum value that progressively decreases with time ( Figure 3C). This pattern is a consequence of binary division and the assumption that the critical C content when cells divide, C max , is equal among cells.
Together, these results show that cell division causes the timedependence of x E S in a cell modeled by zero-order kinetics to follow a zig-zag pattern that closely follows the exponential function in Eq. 15 describing the time-dependence of x E S in a cell modeled by first-order kinetics (Figures 2A, 3A,B, 4A). In contrast, the assumption of zero-order kinetics of C assimilation with no accounting for cell division leads to a dramatically different pattern, where x E S approaches 1 at a much slower rate and the variability among cells, caused for example by intercellular heterogeneity in the C content at the time of the labeled substrate addition, remains pronounced over longer time scales (compare dashed and solid lines in Figure 3A).

Calculating Rates of Cellular C Assimilation
Using the findings presented above, we suggest a three-step procedure for calculating the cellular rates of C assimilation based on 13 C atom fractions measured by nanoSIMS (Figure 5). We emphasize the underlying assumptions of each step to provide the foundations of the formulas. Numerical implementation of the calculation steps, including the calculation of the best estimate and uncertainty of the rates, is provided as a Matlab script that can be interfaced with experimental data organized in a spreadsheet. The script is available at https://github.com/ lpolerecky/LARS. Previous studies documented that isotopic composition of cells can be significantly affected by label loss or dilution associated with common techniques of sample preparation for nanoSIMS analysis including chemical fixation, dehydration and resin embedding (Musat et al., 2012;Hoppe et al., 2013;Loussert-Fonta et al., 2020). If needed, these effects can be corrected for as previously described (Stryhanyuk et al., 2018) before applying the steps described below to calculate assimilation rates.
Step 1 converts the measured 13 C atom fraction of the cell, x, to the source-normalized excess 13 C atom fraction, x E S (Eq. 10). This step accounts for variability in 13 C labeling of the cell due to the 13 C-labeled source, x S,eff (see Eq. 6). Quality results require that x S,eff be well constrained by the experimental set-up, such as by direct measurement. The initial 13 C atom fraction of the cell, x i , is obtained by nanoSIMS measurements of control (unlabeled) cells that were prepared in the same manner as treated (labeled) cells.
Step 2 calculates the carbon-specific rate of C assimilation, k, according to where x E S is obtained in Step 1. This calculation accounts for variability in 13 C labeling of the cell due to the incubation time (Eq. 15). The underlying assumption of this step is that the measured cell assimilated C during the SIP incubation according S in a population of cells with partially synchronized cell cycles. All cells assimilated C at the same cell-specific rate, r. Simulations were performed for the same values of C and r as those shown in Figure 2. The distribution of C i among cells is described by Eq. 18 with C 0 = 9 fmol C cell −1 (red) or C 0 = 13 fmol C cell −1 (blue) and the degree of synchronicity γ = 1.73 (corresponding to C = 1 fmol C cell −1 , both cases). to first-order kinetics or, equivalently, that the measured cell represents the average cell in a population that assimilated C according to zero-order kinetics, had the same cell-specific rate r, and had perfectly unsynchronized cell cycles ( Figure 3B). Importantly, this step accounts for the possibility that the measured cell is the product of cell division that occurred during the SIP incubation.
Step 3 calculates the cell-specific rate of C assimilation, r. Because r is an absolute measure of the C assimilation rate (mol C cell −1 h −1 ), its accuracy depends on the knowledge of the C content of the measured cell (Supplementary Methods, Section "Estimating Cellular C Content"). There are three approaches to calculate r. Approach A applies when the C content of the measured cell is not known precisely but can be approximated by the average C content of the cell over a cell cycle, C (Figure 5, Step 3A). In this case, r is calculated according to using the value of k obtained in Step 2. Approach B applies when only the C content of the measured cell, C, but not the average C content C , can be constrained (Figure 5, Step 3B). In this case, r is calculated according to using the value of k obtained in Step 2. Finally, Approach C applies when both C and C can be constrained (Figure 5, Step 3C). In this case, r is calculated according to using the value of x E S obtained in Step 1. Here, Z −1 is the inverse of the zig-zag function introduced above (Figure 2A), and C max is calculated from C using Eq. 17. We emphasize that all three approaches account for the possibility that the measured cell is the product of cell division that occurred during the SIP incubation.
The assumption underlying Approaches A and B is that the measured cell assimilated C according to first-order kinetics (i.e., at a constant k) during the SIP incubation. Because C in Eq. 21 is the average C content of the measured cell across the cell cycle, r calculated by Approach A represents the average cellspecific rate across the cell cycle. Thus, inherent to Approach A is the assumption that the measured cell assimilated C according to first-order kinetics across the entire cell cycle and with the same value of k as determined from the SIP experiment (i.e., by Eq. 20). In contrast, r calculated by Approach B represents the instantaneous cell-specific C assimilation rate at the time of sampling. Note that the assumption of first-order kinetics implies that r would have been changing during the SIP experiment (because r = k · C, where k is constant, and C is time-dependent; Figure 2E). Thus, r calculated by Approach B represents a rate at a specific time point during the cell cycle of the measured cell, namely at the end of the SIP experiment, when the C content reached the value of C used in Eq. 22.
The assumption underlying Approach C is that the measured cell assimilated C according to zero-order kinetics (i.e., at a constant r) during the SIP incubation and that the C content across the cell cycle varies strictly between the critical values of C max /2 and C max . Because this approach accounts for the possibility that the measured cell is the product of cell division that occurred during the SIP incubation, r calculated by Approach C represents the average cell-specific rate during the SIP incubation. Note that because Approach C makes use of both C and C , it allows the reconstruction of the cell cycle stage of the measured cell. Thus, if the SIP experiment is conducted over a time interval that is short compared with the doubling time of cells in a population, approach C can potentially reveal how the cell-specific rates vary across the cell cycle.

Application to Hypothetical Data
To illustrate the utility of the procedure described above, we calculated C assimilation rates based on hypothetical data from FIGURE 5 | Flowchart for quantifying C assimilation rates in single cells. First, cells are incubated with a 13 C-labeled substrate and their 13 C atom fractions, x, are measured by nanoSIMS. Second, the C content of the cells is constrained. Rate calculation then proceeds in three steps: (1) convert x to x E S to account for the isotopic composition of the substrate, (2) calculate the carbon-specific C assimilation rate, k (this rate is independent of the C content of the cell), and (3) calculate the cell-specific C assimilation rate, r. Depending on how well the C content of the measured cell can be constrained, the third step proceeds either by using the value of k from step 2 (Approach A and B) or using the value of x E S from step 1 and the inverse of the zig-zag function, Z −1 (Approach C). The calculated r then represents the average rate across the cell cycle (Approach A), the instantaneous rate at the end of the SIP experiment (Approach B), or the average rate during the SIP experiment (Approach C). Steps 2 and 3 consider that the measured cell may be a product of a cell that divided during the incubation. The calculation approaches are implemented in a Matlab script available via GitHub (https://github.com/lpolerecky/LARS). cells incubated with a 13 C-labeled C source (x S,eff = 0.1) for 2 hours. Measured cells had low, intermediate, and high levels of 13 C labeling (x 1 = 0.02, x 2 = 0.033, x 3 = 0.09). Using Steps 1 and 2, these values yield carbon-specific assimilation rates of k 1 = 0.053 h −1 , k 2 = 0.14 h −1 , and k 3 = 1.1 h −1 , which correspond to doubling times of about τ 1 = 13 h, τ 2 = 4.9 h, and τ 3 = 0.63 h, respectively.
To calculate the cell-specific rate, r, we further assumed three values for the C content of the measured cell (C = 7.6, 10.4 and 13.2 fmol), each quantified with the same analytical precision of C = 0.1 fmol, and an average C content of C = 10 fmol (which corresponds to the critical C content of C max = 13.9 fmol; Eq. 17). Figure 6 shows that r values predicted using Approach A do not depend on the C content of the measured cell (black bars), whereas they increase proportionally to the C content when calculated using Approach B (gray bars). This pattern is expected from Eqs. 21 and 22. In contrast, r calculated by Approach C (open bars) does not necessarily vary monotonously with the C content of the measured cell. This occurs if the number of divisions during the SIP incubation differs between cells with the same 13 C labeling but with different C contents at the end of the incubation. For example, the cell with C = 7.6 fmol divided once, whereas the cell with C = 10.4 fmol did not divide during the SIP incubation ( Figure 6B, open bars). Thus, the first cell needs to assimilate C more rapidly than the second cell so that both reach the same 13 C labeling at the end of the incubation. Differences between the r values calculated by Approach A and C tend to diminish with increasing 13 C labeling of the cells (compare black and open bars in Figures 6A-C). This pattern emerges because the difference between the first-order and zero-order models of C assimilation decreases with the number of cell divisions during the incubation and thus with the increasing extent of 13 C labeling (Figure 3).

DISCUSSION
This work uses mathematical modeling to identify and explicitly formulate fundamental assumptions underpinning calculations of substrate assimilation rates into cells based on SIP-nanoSIMS data. The 3-step approach to calculate substrate assimilation rates presented in this research evaluates the impact of the model of C assimilation by a cell (zero-order or first-order kinetics), cell division, and the carbon content of the cell measured by nanoSIMS on the calculated rate. Below, we summarize the key assumptions underlying this 3-step approach and compare it to previously published approaches for determining rates of substrate assimilation. Then, we discuss how the rates evaluated in single cells can be upscaled toward a cell population. Finally, motivated by results from a recent nanoSIMS study on diazotrophic cyanobacteria (Polerecky et al., 2021), we identify key sources for the variation in calculated assimilation rates typically observed among cells of the same species within a population that have not been previously discussed.

Summary of Assumptions Used to Calculate Assimilation Rates
Important assumptions underlying the calculations described in this study are: (1) Catabolism or turnover of cell material, such as protein and carbohydrates, does not occur even though their half-lives may be shorter than the incubation duration. Thus, the cell-specific C assimilation rate directly reflects the increase in the cell biomass over time (Eq. 1).
(2) The cell's activity is constant during the incubation. That is, the C assimilation rates r and k in the zero-order and first-order kinetic models, respectively, are time-independent. (3) Carbon density of the cells is constant or unaffected by accumulation of storage inclusions. (4) A cell divides into two identical cells when its C content reaches a certain critical value C max . (5) The effect of kinetic isotope fractionation is minor and negligible. (6) Isotopic composition of the substrate, x S,eff , does not change during the incubation. Analysis of a SIP incubation with timedependent x S,eff is discussed in more detail in Supplementary Methods, Section "Modeling Cellular Assimilation of Substrates With Time-Dependent Isotopic Labeling."

Comparison With Previous Studies
Commonly used approaches for calculating substrate assimilation rates (Foster et al., 2013;Krupke et al., 2015;Svedén et al., 2015;Klawonn et al., 2016;Stryhanyuk et al., 2018;Calabrese et al., 2019;Mills et al., 2020) assume a zero-order kinetic model of substrate assimilation and do not consider that the measured cell could be a product of cell division during the SIP incubation. Additionally, they typically use an average value for the C content or the C density combined with the biovolume of the measured cell. Here we evaluate the impact of these assumptions on calculated rates.
The most common published approaches calculate r according to formulas equivalent to (Krupke et al., 2015;Svedén et al., 2015;Klawonn et al., 2016;Mills et al., 2020;Ploug, 2021) or (Stryhanyuk et al., 2018;Calabrese et al., 2019) In these formulas, C i and C f denote the C content of the measured cell at the beginning and end of the SIP incubation, respectively. Equations 24 and 25 are equivalent and follow directly from the mass balance of a non-dividing cell. Equation 24 is obtained by defining the cell-specific C assimilation rate as r = C a /t, where C a is the amount of carbon assimilated by the cell during time t, and by considering that x E S = C a /C f for a cell that did not divide during time t (see Supplementary Methods, Section "Relationships between x E S , X net , Fx net and K A "). Similarly, Eq. 25 derives from the same assumptions and by additionally considering that, for a cell that did not divide, C i and C f are related according to C f = C i + C a . In summary, these equations yield the average cell-specific C assimilation rate during the SIP incubation. However, the calculated r is accurate only if the cell did not divide during the incubation and if C f and C i are stringently constrained.
Constraining C i or C f of a cell is not trivial even without considering cell division. Facing the practical difficulties of their direct measurement, some practitioners use C in place of C f (e.g., Svedén et al., 2015;Olofsson et al., 2019a) or approximate C i using C f (Stryhanyuk et al., 2018;Calabrese et al., 2019). These substitutions, however, invalidate the approach based on Eq. 24 or 25 because the mass balance underpinning these formulas is no longer achieved.
Many studies estimate the biovolume of the measured cells from nanoSIMS images and apply literature values for the C density, or an empirical relationship between the C density and biovolume of phytoplankton cells (e.g., Verity et al., 1992;Stryhanyuk et al., 2018;Khachikyan et al., 2019), to approximate C f (e.g., Foster et al., 2013;Krupke et al., 2015;Schoffelen et al., 2018;Mills et al., 2020;Trembath-Reichert et al., 2021). Although this approach may yield a well-constrained value of C f , one must be careful when using it in Eq. 24 to calculate r. Specifically, if the cell divided during the SIP incubation, the correct value of C f should include the C content of the measured cell as well as its sister cell (or sister cells, if the cell divided more than once). Only then will the mass-balance requirement be fulfilled, and Eq. 24 will yield the correct value of r for the measured cell. This could, in principle, be achieved by including a correction factor when calculating the cell-specific rate using Eq. 24. Several recent studies (Schoffelen et al., 2018(Schoffelen et al., , 2019Olofsson et al., 2019a,b;Ploug, 2021) have employed a factor 2 in the calculation of carbon-specific growth rates. However, as this latter quantity is solely determined by the cell's isotopic composition (i.e., it does not describe absolute amounts of carbon assimilated; see Eq. 20), FIGURE 6 | Cell-specific rates of C assimilation in cells determined using hypothetical data. The rates were calculated using the procedure described in the text (Steps 1-3) and using approaches A, B and C for Step 3. Each panel (A-C) corresponds to a different value of the 13 C atom fraction of the cell, x, as shown above the graph. The rates represent the average rate across the cell cycle (Approach A), the instantaneous rate at the end of the SIP incubation (Approach B), and the average rate during the SIP incubation (Approach C). Values were calculated for C = 10 fmol C cell −1 , x S,eff = 0.1, t = 2 h, and three values of the C content of the measured cell (x-axis), each determined with the precision of 0.1 fmol. Numbers in brackets indicate how many times the measured cell divided during the incubation, as determined by the algorithm implementing Approach C. Note differences in the y-scale. this operation does not correct for the carbon content of a sister cell as described above.
To illustrate quantitatively the impact of cell division on calculated r, we compared the results obtained using Eq. 24 and our Approach C (Eq. 23) for a range of 13 C labeling and C content of the measured cell (Figure 7). We did not apply a correction factor to calculate C f but assumed that C f in Eq. 24 is the same as C in Eq. 23. As expected, both approaches yield the same result if the cell did not divide during the SIP incubation (see values in Figure 7 for which the number of cell divisions is equal to 0), because the function in Eq. 11, which is the basis for Eq. 24, is equal to the zig-zag function Z (see Figure 2A, black curve before t/τ < 1 and red curve before t/τ < 0.5). However, if the measured cell did divide during the incubation, Eq. 24 yields a lower r value compared to that calculated by Approach C (Figure 7). The difference between the two approaches becomes more pronounced as x E S of the measured cell increases above a threshold value that depends on the cell cycle stage (Figure 7B). For example, the threshold value of x E S is 0.1 if the cell cycle stage of the measured cell is 10%, but it increases to 0.33 and 0.47 for the cell cycle stage of 50% and 90%, respectively (inset in Figure 7B).
The question remains, however, how to determine whether a cell measured by nanoSIMS has divided during the SIP incubation. Tracing the history of a cell analyzed by nanoSIMS is practically impossible because nanoSIMS is a terminal measurement. The history can be reconstructed if additional information about the measured cell, namely C , can be constrained. Constraining C allows evaluation of the minimum and maximum C content of the measured cell (C max /2 and C max , respectively; Eq. 17). Combining this information with the carbon content, C, enables estimation of the cell cycle stage of the measured cell. Combined with the value of x E S , the history of the measured cell during the SIP experiment can be reconstructed to estimate the value of r that accounts for cell division. This sequence of steps is implemented by our Approach C based on the inverse of the zig-zag function (Eq. 23).
We conclude that the 3-step procedure proposed in this study (Figure 5) is applicable to many situations encountered by users of SIP-nanoSIMS to estimate substrate assimilation rates in individual cells. When the C content of a measured cell, C, cannot be constrained, the carbon-specific assimilation rate, k (Eq. 20), and the corresponding doubling time, τ = ln2 /k, is the maximal information that can be estimated from the measured isotopic composition of the cell. Constraining C allows estimation of the instantaneous cell-specific assimilation rate at the end of the SIP experiment (Approach B, Eq. 22), while constraining C allows extrapolation towards the average cell-specific assimilation rate across the cell cycle (Approach A, Eq. 21). Finally, constraining both C and C allows determination of the average cell-specific assimilation rate over the SIP experiment (Approach C, Eq. 23). Note that accounting for cell division, as done in Approach C, critically depends on the assumption that the measured cell divided during the incubation when its C content reached C max estimated from C . If this assumption is not valid (e.g., because the C assimilated by the cell was allocated into storage inclusions, or cell division was delayed or not binary), or cannot be verified, the instantaneous rate at the end of the SIP experiment (Eq. 22) is the maximal information about the cell-specific assimilation rate that can be estimated by combining the 13 C labeling and C content of the measured cell.

Extrapolation Toward a Cell Population
A frequent aim of SIP-nanoSIMS measurements is to estimate a bulk C assimilation rate for a target cell population and thus evaluate its potential impact on C fluxes in the environment. This upscaling can be done by first averaging the cell-specific C assimilation rates determined for individual target cells to estimate the population average, r , and then multiplying r by the total target cell abundance in the environment. When calculating the assimilation rates for individual cells, Approach A or B is recommended if cell cycles in the target population are perfectly unsynchronized, whereas Approach C is recommended for partially synchronized cells. Because these approaches account for cell division, r estimated in this way will reflect the population average across a broad range of incubation times (Supplementary Figure 5), which may be important when using SIP-nanoSIMS to study environmental FIGURE 7 | Comparison of previously published and the current approaches for calculating C assimilation rates in cells. The key difference is that our Approach C accounts for cell division, whereas cell division is not considered in the approach using Eq. 24. (A) Carbon-specific rates predicted from the 13 C labeling of the measured cell, expressed as x E S . Calculations assumed C = 10 fmol C cell −1 , t = 2 h, and three values of the cell's C content, C.
The values of C = 7.6, 10.4, and 13.2 fmol correspond to the cell cycle stage of s = 0.1, 0.5, and 0.9, respectively. (B) Ratio between the r value predicted by Approach C and the approach based on Eq. 24, calculated for several values describing the cell cycle stage of the measured cell (indicated by numbers next to the lines). Approach based on Eq. 24 yields lower rates compared to Approach C if the measured cell divided during the incubation [see number of divisions in the top graph in panel (A)]. This occurs when x E S exceeds a threshold value that increases with the increasing cell cycle stage of the measured cell (see inset in panel (B)). Note that the ratio only depends on x E S and s and not on the incubation time.
samples where the growth rates of cells are a-priori unknown.
Averaging of cell-specific rates calculated using Eq. 24 is not recommended if the expectation is that a significant number of target cells divided during the incubation, because this approach would underestimate r (Supplementary Methods, Section "Simulating SIP Incubation, " Supplementary Figure 5). Note that the precision and accuracy of the estimated bulk assimilation rate will still depend on the number of target cells measured by nanoSIMS, i.e., on the quality of sampling of the target population.

Interpretation of Calculated Rates and Their Variability Among Cells
Typically, 13 C atom fractions obtained by nanoSIMS measurements vary among cells, raising the question, to what extent this variability is caused by differences in the intrinsic metabolic activities of cells probed by the SIP experiment (i.e., r) versus differences arising from cells at different stages of their cell cycle during the labeling interval. Additionally, we note that assimilation rates obtained in nanoSIMS studies can reflect variation not only in the targeted process (e.g., C fixation) but also in other metabolic processes and cellular characteristics (e.g., simultaneous assimilation of other external or internal C sources). In this study we identified that in addition to r, the sources of variability include t, x S,eff , and C i . Our 3-step procedure ( Figure 5) accounts for some, but not all, of these sources of variability. Below, we discuss what this uncertainty implies for interpreting the calculated rates of C assimilation (k or r) and especially their variability among measured cells. We focus on C, but similar arguments can be made for other elements.

Influence of an Alternative C Source
The effective isotopic composition of the C source, x S,eff , will almost certainly be different from the isotopic composition of the target C source, x S , if cells assimilate carbon from an alternative source in addition to the target source. The parameter f tar accounts for the influence of alternative C sources (see Materials and Methods and Table 1). Calculated k values will underestimate total C assimilation if an alternative C source influenced x S,eff ( Table 2). Thus, cell-to-cell differences in k may reflect intercellular variation in use of the target source rather than total C assimilation. Cell-specific rates, r, will be similarly affected as they are related to k.

Influence of C Storage Content
Because k is defined as k = r/C, variability in calculated k among cells can be caused not only by variation in the instantaneous C assimilation rate (r) but also by variation in the instantaneous C content (C). Such variations in C can be expected, for example, in cells that encounter starvation, i.e., they have depleted a required growth resource from the environment, or the environment has Shown are results for two cells that assimilate C with equal total carbon-specific rates of 0.01 h −1 over 2 h. The cells differ only in the fraction of C assimilated from the target C source, f tar . If C is incorrectly assumed to be assimilated from only the target source, calculated k will underestimate the total C assimilation rate (as discussed in Section Influence of an Alternative C Source). The calculated and total values will only be equal if C assimilation from the alternative C source is accounted for by using the correct value for x S,eff (lines 3-4, values shown in bold). Calculated for x S,tar = 0.99, x S,alt = 0.011, and x i = 0.011. changed, such that assimilation of that resource is no longer possible. These conditions cause cells to cease exponential growth and commonly lead to the synthesis of storage compounds (e.g., polysaccharides). The variability in calculated k (Eq. 20) among starving cells may therefore reflect variation in storage content rather than r.
Similarly, cell-to-cell variability in calculated r needs to be interpreted with caution as it may reflect variation in storage content, rather than intrinsic cell-specific C assimilation rates, because neither Approach A, B nor C accounts for potential variation in the cell's C content caused by storage (see the underlying assumptions in Figure 5). Specifically, r is calculated either directly from k using C (Approach A, Eq. 21), or the calculation relies on the assumption that the cell C density, which may be used to estimate the C content of the measured cell based on its biovolume (Supplementary Methods, Section "Estimating Cellular C Content"), is constant during the SIP experiment and equal among cells (Approach B and C,. These assumptions are unlikely to be valid in cells that metabolize storage inclusions. By accounting for the relationships between storage content and apparent assimilation rates, however, we can deduce information about storage metabolism from nanoSIMS data, as outlined in the following section.

Interpretation of Simultaneous 13 C and 15 N Labeling Experiments
For a single labeled element (e.g., 13 C or 15 N), the importance of storage metabolism versus assimilation of new material is difficult to resolve without additional measurements at subcellular level (e.g., TEM images of storage inclusions). On the other hand, dual labeling (i.e., simultaneous incubations with two labeled substrates such as 13 C and 15 N) can greatly increase informative value of the SIP experiment. By systematic consideration of the effects of storage inclusion metabolism on the element-specific rates of C and N assimilation (k C and k N ), whole-cell nanoSIMS measurements can be used to gain information not only about the acquisition of new carbon and nitrogen, but also about intracellular fluxes of these elements (Polerecky et al., 2021). Specifically, metabolic pathways leading to storage inclusion biosynthesis and catabolism can be revealed by leveraging information about the known elemental compositions of specific storage inclusions (e.g., polysaccharides or cyanophycin) and the isotopic compositions of target C and N sources in carefully designed SIP experiments.
We simulated assimilation of 13 C-and 15 N-labeled substrates by cells under various scenarios, including preferential incorporation of the assimilated C and N into different cell compartments (e.g., cell matrix or storage inclusions), or C and N incorporation using internally recycled C and N in addition to C and N originating from the target sources provided externally (Supplementary Methods, Section "Modeling Simultaneous Assimilation of C and N"). Subsequently, we calculated the Cand N-specific assimilation rates, k C and k N , in these cells using the approach described above (Eq. 20).
If a cell assimilated the 13 C-and 15 N-labeled source substrates in a balanced way, i.e., the amounts of assimilated C and N did not change the C:N ratio of the cell, the ratio of the calculated element-specific rates of C and N assimilation, k C /k N , is equal FIGURE 8 | Simulations of cell-specific assimilation of C and N. Element-specific rates of C and N assimilation by whole cells, k C and k N , calculated from hypothetical data. The first simulation (stars) shows the results for assimilated C and N that originated from two sources: labeled target sources [x( 13 C) S,tar = 1, x( 15 N) S,tar = 1] and unlabeled alternative sources [x( 13 C) S,alt = 0.011, x( 15 N) S,alt = 0.0037]. Resultant k C and k N values assumed assimilation of only the target sources [x( 13 C) S,eff = 1, x( 15 N) S,eff = 1]. Individual data-points correspond to cells that assimilated different fractions of C and N from the target sources, with corresponding f C,tar and f N,tar values given in parentheses (as f C,tar :f N,tar ). Deviation of k C /k N from 1 reveals assimilation of C or N from alternative sources. The second simulation shows the impacts of preferential incorporation of C and N into cyanophycin granules (cy, circle), polysaccharide inclusions (ps, triangle), and cell matrix (m, square). Cellular compartments with C:N that differ from the whole cell can cause k C /k N to deviate from 1. The third simulation (diamonds) shows results for cyanophycin synthesis using unlabeled C and N from inside the cell (e.g., protein and polysaccharide catabolism) in addition to C and N from the external target sources. The parentheses give f C,tar :f N,tar . Use of unlabeled internal C and N pools causes the cell to assimilate C and N from sources with x( 13 C) S,eff < 1 and x( 15 N) S,eff < 1, whereas calculation of k C and k N assumed x( 13 C) S,eff = 1 and x( 15 N) S,eff = 1. Details of the simulations are given in Supplementary Methods, Section "Modeling Simultaneous Assimilation of C and N." to 1. Deviation of k C /k N from 1 can arise for multiple reasons. First, the effective isotopic composition of the C or N source, x( 13 C) S,eff or x( 15 N) S,eff , differed from the isotopic composition of the target C and N source (Figure 8, stars). For example, by providing autotrophs with 13 C-DIC (their only carbon source) and 15 N-ammonium, a value of k C /k N > 1 could indicate uptake of organic N compounds in addition to ammonium. Thus, deviation of k C /k N from 1 can reveal uptake of additional substrates. A similar approach was used by Schoffelen et al. (2018) to reveal uptake of organic P in addition to inorganic P. Second, C and N assimilation did not match the C:N ratio of the whole cell. For example, 13 C or 15 N were incorporated into inclusions with C:N ratios that differed from that of the whole cell (Figure 8, circle, square, triangle). This can happen in cells that temporally decouple assimilation of C and N, such as in some diazotrophs, which assimilate C in the light and N mostly in the dark (e.g., Polerecky et al., 2021). Third, C or N from existing cell material (e.g., proteins, storage compounds) was mobilized into a new form. For example, protein degradation released N that was then used for cyanophycin synthesis in addition to incorporation of 15 N from the target source (Figure 8, diamonds; see also Polerecky et al., 2021). Thus, deviation of k C /k N from 1 can be used to identify the presence of intra-cellular C and N stores and evaluate their synthesis or mobilization without the need for direct measurement at the sub-cellular level. Notably, the assessment of the deviation of k C /k N from 1 is solely based on the nanoSIMS measurement and does not require knowledge about the C and N content of the measured cell.

CONCLUSION
This study revisited calculations of biologically relevant parameters from SIP-nanoSIMS data. We suggest a step-bystep procedure to evaluate elemental assimilation rates into whole cells, discuss key assumptions underlying the procedure, and describe factors that can introduce variability into the calculated rates.
Determining single-cell assimilation rates is challenging, as it requires proper evaluation of multiple factors that are difficult to constrain experimentally or analytically, necessitating assumptions in the calculation process. Our analysis illustrates and quantifies the consequences of some of the key assumptions, particularly focusing on assimilation kinetics (zero vs. first order), cell division, and the elemental content of the measured cell. We give guidance for interpreting the rates calculated by approaches that differ depending on knowledge of the elemental content of the measured cell and develop an approach to account for cell division in the calculation process.
We call for caution when interpreting inter-cellular variability of assimilation rates calculated from SIP-nanoSIMS data if (1) cells can assimilate elements such as C from alternative sources present in the environment in addition to the target source added externally during the SIP experiment, (2) storage inclusions are synthesized, catabolized or merely present inside the studied cells during the SIP experiment, or (3) cell divisions during the SIP experiment may be expected. Our modeling analysis further shows that for dual labeling SIP experiments (e.g., using 13 C and 15 N labeled substrates), deviation of k C /k N from 1 can be used to identify the presence of intracellular C and N stores and evaluate their synthesis or mobilization without the need for direct measurement at the sub-cellular level.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.