DNA stable isotope probing employs isopycnic centrifugation to separate isotopically enriched (“heavy”) DNA molecules from unlabeled (“light”) DNA based on their differences in buoyant density (BD). Isopycnic centrifugation is distinguished from other centrifugation methods in that centrifugation is carried out long enough to both generate a density gradient (typically using CsCl for DNA-SIP) and allow all macromolecules of interest reach sedimentation equilibrium, which is the point at which sedimentation rates equal rates of diffusion (Hearst and Schmid, 1973; Birnie and Rickwood, 1978). Empirical studies have shown that the average BD (ρ) of a mixture of DNA molecules is linearly related to the average G + C content for that collection of molecules:

where [G + C] is the mole fraction of genome G + C content (Schildkraut et al., 1962; Birnie and Rickwood, 1978). In addition, empirical studies have also shown that homogeneous mixtures of DNA molecules form a Gaussian distribution in an isopycnic gradient when at sedimentation equilibrium (Meselson et al., 1957; Fritsch, 1975). Therefore, in order to model the BD distribution of a heterogeneous set of genomic DNA fragments, a Gaussian distribution must be estimated for each homogeneous subset of molecules rather than using discrete BD values (as described in Supplementary Material). Based on the work of Meselson et al. (1957), Fritsch (1975) derived an equation describing time to reach sedimentation equilibrium, which can be reworked to calculate the standard deviation (σ) of the Gaussian distribution:

where L is the effective length of the gradient (cm), t is time in seconds, ω is the angular velocity (radians sec⁻¹), r_p is the distance of the particle from the axis of rotation (cm), s is the sedimentation coefficient of the particle, β° is the coefficient specific to the density gradient medium (e.g., CsCl); p_p and p_m are the maximum and minimum distances between the gradient and axis of rotation (cm) (Fritsch, 1975). By assuming that sedimentation equilibrium has been reached for all macromolecules of interest, Clay and colleagues derived a simplified equation for determining σ from the calculations in Schmid and Hearst (1972):

where ρ is the BD of the particle, R is the universal gas constant, T is the temperature in Kelvins, β is a proportionality constant for aqueous salts of specific densities, G is a buoyancy factor as described in Clay et al. (2003), M_C is the molecular weight per base pair of DNA, and l is the fragment length (bp). For most DNA-SIP experiments, the assumption of sedimentation equilibrium for all DNA fragments is likely to be unrealistic for relatively short DNA fragments (e.g., <4 kb), given that the time to reach equilibrium is inversely proportional to diffusion and hence rises dramatically with decreasing fragment length (Meselson et al., 1957; Birnie and Rickwood, 1978; Youngblut and Buckley, 2014). However, the ultracentrifugation durations used in typical DNA-SIP experiments should still generally produce small σ values for short DNA fragments according to Equation (2) (Neufeld et al., 2007). Therefore, Equation (3) provides a good approximation for modeling the BD distribution of DNA in density gradients generated in typical DNA-SIP experiments.

The distribution of a heterogeneous mixture of DNA fragments in an isopycnic gradient can thus be modeled by integrating the Gaussian distributions of each homogeneous subset of DNA fragments, where the mean of each Gaussian is determined by Equation (1) and the standard deviation derived from Equation (3). In this way, the BD distribution for a given genome in an isopycnic gradient can be modeled by the following steps: simulate genome fragmentation resulting from DNA extraction, bin gDNA fragments with respect to length and G + C content, model Gaussian distribution for each fragment bin, and then integrate these distributions to describe the cumulative DNA distribution in fractions of the gradient. A strictly Gaussian model predicts that, within a CsCl density gradient, DNA fragments should be undetectable (i.e., probability density < 10⁻⁷) in factions of both high and low BD (Figure S1). However, empirical observations show DNA to occur throughout CsCl gradients (Birnie and Rickwood, 1978; Lueders et al., 2004; Leigh et al., 2007). We were able to reconcile the difference between observed and expected DNA distributions as a function of fluid mechanics during gradient reorientation (as described below and in Supplementary Material).

Gradient reorientation impacts the BD distribution of DNA. During isopycnic centrifugation, the buoyant density gradient forms perpendicular to the axis of rotation (Figure S2), and gradient reorientation during centrifuge deceleration is dramatic, especially for vertical rotors (Flamm et al., 1966). While the distortion of the BD gradient during reorientation has been shown to be minimal in the aggregate (Fisher et al., 1964; Flamm et al., 1966), the inevitable presence of a diffusive boundary layer along the tube wall is sufficient to entrain quantities of DNA, which are small but should be readily detectable by high throughput sequencing methods. The flow field that occurs during gradient reorientation entrains along the tube wall a volume with a dimension proportional to flow velocity, fluid viscosity, and surface topography (Tritton, 1977; Cohen and Dowling, 2012). Following gradient reorientation, DNA from the entrained volume will combine with DNA from the reoriented volume, thereby introducing a small amount of non-BD-equilibrium DNA into each gradient fraction (Figure S2). The ability of the diffusive boundary to introduce non-BD-equilibrium DNA into gradient fractions can be modeled as a function of rotor geometry and the effect is more pronounced in vertical rotors relative to fixed angle rotors (Figure S2). Assuming sedimentation equilibrium, BD (ρ) can be directly related to the distance from axis of rotation (Birnie and Rickwood, 1978):

From this calculation, the location of DNA molecules in the centrifuge tube, both during centrifugation and fractionation, can be ascertained by using simple trigonometry along with knowledge of centrifuge tube dimensions and angle to the axis of rotation. A full description of the calculations along with an example can be found at https://github.com/nick-youngblut/SIPSim. The fraction of a taxon's DNA fragments that are in the boundary layer (D_ti) is modeled as:

where A_ti is the pre-fractionation community relative abundance of taxon t in gradient i, γ is a weight parameter determining the contribution of A_ti to A_b, and α is the baseline fraction DNA in A_b.

Assimilation of the commonly used isotopes ¹³C and ¹⁵N into genomic DNA produces linear shifts in BD, with a maximum shift of 0.036 and 0.016 g ml⁻¹, respectively (Birnie and Rickwood, 1978). Thus the shift in BD (ρ) can be modeled as:

where I_i,max is the maximum possible BD shift if 100% atom excess for isotope i, A is the atom % excess of isotope i, and ρ_12C is the buoyant density at 0% atom excess. Note that the same formula applies to ¹⁵N.

Based on the theory described above, our SIP data simulation framework simulates the distribution of gDNA fragments in isopycnic gradients at sedimentation equilibrium. Furthermore, it generates the DNA-SIP datasets obtained from fractionating isopycnic gradient(s) and performing high throughput sequencing on many of the gradient fractions. Our framework also implements all of the DNA-SIP analysis methods assessed in this study (Heavy-SIP, HR-SIP, MW-HR-SIP (see below), qSIP, and ΔBD) and evaluates their accuracy of identifying incorporators or quantifying BD shifts. An overview of our simulation framework is shown in Figure 1.

Our simulation framework is a modular collection of steps that can be grouped in workflow stages that are further broken down into steps (Figure 1). The input is a set of reference genomes in fasta format and a text file designating the experimental design, which includes the number of gradients for labeled treatments and unlabeled controls.

Stage 1 involves generating a BD distribution of gDNA fragments for each genome. Step 1a involves simulating the pool of gDNA fragments that is extracted from SIP incubation samples and then loaded into the isopycnic gradients. If amplicon sequence data (e.g., SSU rRNA) is to be generated, amplicons from only the fragments containing the PCR template (“amplicon-fragments”) are sequenced, while shotgun metagenomic sequencing can target all gDNA fragments (“shotgun-fragments”). If ≥1 PCR primer set is provided, amplicon-fragments are generated from genomic regions fully encompassing genome locations that produced amplicons by in silico PCR. Alternatively, shotgun-fragments are randomly generated from all possible genomic locations. The fragment size distribution is user-defined (Table S1).

As described in Equations (1, 3), the length and G + C content of a DNA fragment can be used to calculate a probability distribution of its location in the gradient, assuming sedimentation equilibrium. Step 1b uses the fragments simulated in Step 1a to generate a 2-dimensional Gaussian kernel density estimation (KDE) for each taxon, which describes the joint probability of obtaining fragments with a certain length and G + C content from that taxon. From this 2D-KDE, a large number of [length, G + C] vectors can be simulated efficiently for more precise estimations of the fragment BD distributions. Fragment BD distributions are calculated for each taxon in Step 1c by sampling [length, G + C] vectors from the 2D-KDE and calculating Gaussian distribution from each, where the mean is based on Eq. 1 and the standard deviation based on Equation (3). The collection of Gaussian distributions for all fragments for each taxon is integrated into a BD distribution for all fragments of a taxon with Monte Carlo error estimation, which involves sampling BD values from the collection of Gaussian distributions and estimating a probability density function (PDF) of the fragment BD distribution as a one-dimensional Gaussian KDE. The result is a list of KDEs, with each describing the probability of detecting the gDNA fragments of a taxon at any point along the isopycnic gradient. These fragment BD distributions are modified in steps 1d and 1e by adding diffusive boundary layer (DBL) effects (see section Theory Underlying the Simulation Framework) and isotope incorporation, respectively. The “smearing” due to DBL effects is modeled as a uniform distribution describing the increased fragment BD uncertainty, and this uncertainty is integrated into the fragment BD distributions by Monte Carlo error estimation as in Step 3b. The BD shift due to isotope incorporation is modeled in a similar manner, except BD uncertainty is a result of inter- and intra-population variation in the amount of isotope incorporated. Variation of isotope incorporation is modeled as a hierarchical set of mixture models (weighted sets of standard distributions; such as two Gaussians), where the parameters for intra-population mixture models that describe the amount of isotope incorporated by each individual are themselves defined by inter-population mixture models that describe how isotope incorporation varies among taxa.

Stage 2 involves simulating the isopycnic gradients for a particular experimental design. Step 2a involves simulating the BD range size of each fraction of each gradient. Sizes are drawn from a user-defined distribution. Step 2b involves simulating the relative abundance distribution of taxa in the gDNA pools loaded into each gradient (“pre-fractionation communities”). The abundance distribution of each pre-fractionation community is user-defined and can vary among gradients. Furthermore, the amount of taxa shared or rank-abundances permutated among communities (i.e., the beta-diversity) is user-defined.

Stage 3 involves generating a DNA-SIP dataset based on the fragment BD distributions simulated in Stage 1 along with the isopycnic gradient data generated in Stage 2. In Step 3a, an OTU (taxon) abundance table is generated by sampling from the fragment BD distributions of each taxon generated in Stage 1, with sampling depth determined by pre-fractionation community abundances simulated in Step 2b. The subsampled fragments are then binned into gradient fractions simulated in Step 2a. The resulting OTU table lists the number of gDNA fragments of each taxon in each gradient fraction in each gradient. If the simulated fragments are amplicons, then PCR amplification efficiency biases are simulated in Step 3b based on the PCR kinetic model described in Suzuki and Giovannoni (1996). The model assumes that efficiencies decrease as the product concentration increases due to an increased propensity of single stranded products to re-anneal to their homologous complements. Sequence data is simulated in Step 3c by subsampling from the table of fragment counts (the DNA fragment pool), which produces a final table (“DNA-SIP dataset”) of taxon relative abundances in each gradient fraction in each gradient.

Unless stated otherwise, we made the following assumptions for all simulations in this study. Community abundance distributions were simulated as lognormal distributions with a mean of 10 and a standard deviation of 2. All taxa were shared among communities, and no rank-abundances were permuted. The total number of fragments in each gradient was 1e⁹. Gradient fragment BD range sizes were sampled from a normal distribution, with a mean of 0.004 and a standard deviation of 0.0015. SSU rRNA amplicon-fragments were simulated using the V4-targeting 16S rRNA primers: 515F and 927R (5′-GTGYCAGCMGCMGCGGTRA-3′; 5′-CCGYC AATTYMTTTRAGTTT-3′), as used by Pepe-Ranney et al. (2016a). The amplicon-fragment size distribution was a left-skewed normal distribution with a mean of ~12 kb, which is similar to size distributions produced from common bead beating cell lysis methods (Kauffmann et al., 2004; Roh et al., 2006; Thakuria et al., 2008). A total of 1e⁴ amplicon-fragments were simulated per genome, which equated to >100X coverage for the genomic region of interest. Monte Carlo error estimation was conducted with 1e⁵ sampling replicates. Ultracentrifugation conditions were set as in Pepe-Ranney et al. (2016a), with a Beckman TLA-110 rotor spun at 55,000 rpm for 66 h at 20°C and an average density gradient 1.7 g ml⁻¹. Inter-population variation in isotope incorporation was binary (either 0% or X% atom excess), and intra-population variation was set to zero. Two key parameters were estimated from empirical DNA-SIP data: the bandwidth (smoothing factor) for kernel density estimation, and the gamma parameter in Equation (5). See Table S1 for a full listing of simulation parameters.

The HR-SIP method was performed as described in Pepe-Ranney et al. (2016a,b). Briefly, we used a “heavy” BD window of 1.71–1.75 g ml⁻¹, a sparsity cutoff of 0.25 (i.e., OTUs must be present in >25% of samples), a log₂ fold change null threshold of 0.25, and a false discovery rate cutoff of 10%. ΔBD was determined as described by Pepe-Ranney et al. (2016a), with OTU abundances linearly interpolated across 20 evenly spaced values across the gradient BD range. Briefly, ΔBD is calculated as the difference in the center of mass of abundance distributions (i.e., the average BD weighted by relative abundance) between the labeled treatment and unlabeled control.

qSIP was conducted as described in Hungate et al. (2015), with 90% confidence intervals calculated from 1,000 bootstrap replicates. The variance among qPCR replicates was modeled based on the qPCR data provided in Table S2 of Hungate et al. (2015). Specifically, we found the qPCR count variance (σ²) to increase as a function of the mean (μ). The following polynomial regression was found to best describe this relationship and was used for simulating all qPCR count values:

where μ was set as the total number of simulated DNA fragments in the gradient fraction (designated in the OTU table from Step 4a).

A range of alternative analytical approaches to Heavy-SIP can be used to detect incorporators. We characterized four additional approaches to identifying labeled OTUs. Method 1 identifies as labeled any OTU that occurs in “heavy” fractions of the labeled gradient. Method 2 identifies as labeled any taxa present in the “heavy” fractions of the labeled treatment and absent from the “heavy” fractions of the control gradient. Method 3 identifies as labeled any taxa present in the “heavy” fractions of the labeled treatment and absent in the “light” fractions of the labeled treatment. Method 4 identifies as labeled any taxa present in the “heavy” fractions of the labeled treatment and absent from both the “heavy” fractions of the control and the “light” fractions of the labeled treatment. Of these four approaches, Method 1 provided the highest accuracy (Figure S6) and so this is the method that we used to represent “Heavy-SIP.”

We hypothesized that HR-SIP sensitivity could be improved by altering the “heavy” BD window (1.71–1.75 g ml⁻¹) in which sequence composition is compared between treatment and control. We used SIPSim to evaluate a range of analytical approaches (data not shown) and found that the analysis of multiple windows (hereby called “MW-HR-SIP”) resulted in a significant improvement in sensitivity relative to HR-SIP. MW-HR-SIP evaluates sequence composition within BD windows of: 1.70–1.73, 1.72–1.75, 1.74–1.77 g ml⁻¹ (Figure S3) while adjusting for multiple comparisons.

The genome dataset used to simulate genomic DNA fragments was obtained from Genbank (Benson et al., 2008). From a list of all bacterial genomes designated as “complete,” one representative was chosen per species in order to reduce the bias toward highly represented species. We found the dataset to contain a rather high proportion (~12%) of low G + C organisms (<30% G + C); most of which were obligate endosymbionts. We randomly sampled a subset of these low G + C genomes in order to reduce the proportion of low G + C organisms to just 1% of the genome dataset. The resulting dataset consisted of 1,147 bacterial genomes.

In order to simulate empirical data from Lueders et al. (2004), the genome sequences of Methanosarcina barkeri MS and Methylobacterium extorquens AM1 were downloaded from Genbank. Amplicon-fragments were simulated with the primers Ar109f (5′-ACKGCTCAGTAACACGT-3′), Ar915r (5′-GTGCTCCCCCGCCAATTCCT-3′), Ba519f (5′-CAGCMGCCGCGGTAANWC-3′), and Ba907r (5′-CCGTCAATTCMTTTRAGTT-3′). Atom % excess was assumed to be 100%, and isopycnic centrifugation conditions were simulated as specified in Lueders et al. (2004).

For model evaluation (see Supporting Results), we downloaded the genomes Clostridium ljungdahlii DSM 13528, Escherichia coli 1303, and Streptomyces pratensis ATCC 33331 from Genbank.

We compared the properties of DNA fragment BD distributions between simulated DNA-SIP data and empirical DNA-SIP data obtained from a soil community. The DNA-SIP dataset from Youngblut and colleagues consisted of SSU rRNA MiSeq sequences (V4 region) of ~24 fractions per gradient from 6 gradients of unlabeled controls. This dataset was generated using agricultural soils (15 g per sample) amended with a complex substrate mixture (3.6 mg C per g soil), and incubated aerobically at 50% water holding capacity and at room temperature. DNA was extracted following destructive sampling of replicates at days 1, 3, 6, 14, 30, and 48. DNA was subject to CsCl centrifugation (TLA110 Beckman rotor, 55,000 rpm, 66 h, 1.69 g ml⁻¹ average gradient density) and fractionated (100 μl fractions) using methods which have previously been described in detail (Pepe-Ranney et al., 2016b). Since the use of these samples in the present context is limited to the analysis of DNA buoyant density distribution properties in the gradient, only control samples containing unlabeled DNA were examined. These data were subsampled to obtain a total richness equal to the 1,147 OTUs in our reference genome dataset. The sequence data is available from the NCBI under BioProject PRJNA382302.

The SIP simulation framework was mostly written in Python v2.7.11, with some accompanying code written in C++ v4.9.2 and R v3.2.3 (R Core Team, 2016). MFEprimer v2.0 was used to perform in silico PCR (Qu et al., 2009). The software, along with documentation and examples, can be found at https://github.com/nick-youngblut/SIPSim. All genomes were downloaded from Genbank with the R package genomes v2.12.0 (Stubben, 2014), and all data analysis was conducted in R with the following packages: ggplot2 v2.1.0, dplyr v0.4.3, tidyr v 0.4.1, and cowplot v0.6.2.

Further methodological details are provided in the Supplementary Material.

The SIPSim model starts with a set of user-designated genomes and user-designated experimental parameters (e.g., number of gradient fractions, desired community characteristics, desired isotopic labeling characteristics) as described (see section Methods and Supplementary Material). Briefly, the genomes are fragmented as would occur during DNA extraction, isotopic labeling is applied to some number of genomes as specified by the user, the BD distributions are determined for each DNA fragment and fragment collections are then binned into gradient fractions, fragments are sampled from each fraction as would occur during amplification and DNA sequencing of SSU rRNA genes, and then the relative abundance is calculated for each OTU (Figure 1, see also section Methods). The model produces results that are highly similar to those observed in empirical experiments, including the ability to detect DNA fragments throughout the density gradient (Figure 2).

The development of the simulation model was guided by established centrifugal theory and by comparison of simulated results to empirical data (as in section Methods and Supplementary Material). First, we performed a simple evaluation of model performance by recreating results from a prior DNA-SIP experiment with Methanosarcina barkeri MS and Methylobacterium extorquens AM1 (Lueders et al., 2004) (Figure S4). Simulated DNA distributions (both in terms of total DNA and SSU rRNA gene amplicon copies) significantly and strongly correlated with the empirical data for both taxa (p < 0.003 for all comparisons; see Table S2). In addition, the simulated SSU rRNA gene amplicon-fragment BD distributions were shifted 0.007 g ml⁻¹ toward the middle of the BD gradient relative to the shotgun-fragments (“total DNA”), a phenomenon also observed in the empirical data. This central tendency for SSU rRNA amplicon-fragments reflects G + C conservation of the rrn operon, as previously described (Youngblut and Buckley, 2014).

Next, we evaluated SIPSim results by comparing to empirical data. For this purpose SIPSim output was compared to results obtained with unlabeled DNA from soil (see section Datasets). The empirical data was derived from an experiment in which an unlabeled nutrient mixture was added to soil and DNA was extracted at 1, 3, 6, 14, 30, and 48 days. These six DNA samples were equilibrated in CsCl gradients, fractionated by BD, and SSU rRNA gene amplicons were sequenced for ~24 fractions from each gradient. The simulation run included 1,147 microbial genomes (see section Methods), hence the soil data was resampled to 1,147 OTUs in order to standardize the richness of the simulated and empirical data. The empirical results reveal several interesting features about the distribution of DNA in CsCl gradients. First, we observed that variance in DNA fragment BD is positively correlated with its relative abundance and that most taxa with relative abundances >0.1% are detected in all CsCl gradient fractions (Figure 3). Second, we observed that taxonomic similarity is auto-correlated with fraction BD, so that fractions of similar BD have similar nucleic acid composition (Figure 3). Finally, we observed that changes in community composition caused dramatic shifts in the sequence composition of “heavy fractions” even in the absence of isotopic labeling (Figure 3). We applied three different analytical approaches to assess similarity between empirical and simulated data: correlation in Shannon Diversity with respect to fraction BD, correlation in Jaccard dissimilarity with respect to fraction BD, and correlation in OTU BD range with respect to taxon relative abundance. We found that variance between simulated and empirical results was significantly less than the variance observed between replicate empirical samples (Figure S5). Furthermore, we used these comparisons between simulated and empirical results (Figure S5) to optimize parameter values for use in the SIPSim model (Table S1, as described in Supplementary Material).

We hypothesized that both the number of taxa that incorporate isotope and the atom % excess isotope incorporation per taxon would substantially affect the accuracy of DNA-SIP methods. To test these predictions, we simulated DNA-SIP datasets for both ¹³C-labeled samples and unlabeled controls (3 replicates of each), while varying both the number of incorporators (1, 5, 10, 25, or 50% of taxa) and the atom % excess isotope incorporation for each taxon (0, 15, 25, 50, 75, or 100 atom % excess ¹³C). Taxa in the control were always set to 0% atom excess isotope incorporation. Each simulation was replicated 10 times, with differing taxa randomly designated as incorporators in each replicate. We evaluated 7 methods used to analyze DNA-SIP data including: HR-SIP, qSIP, 4 different approaches to “Heavy-SIP,” and MW-HR-SIP. Heavy-SIP includes a family of approaches (see section Methods) in which incorporators are identified on the basis of presence-absence in heavy and/or light fractions (Figure S6).

The model predicts that both the number of incorporators and the amount of isotope incorporated affected accuracy (Figure 4). However, the predicted effect of these parameters on specificity and sensitivity varied depending on the analytical method (Figure 4). Specificity is the proportion of true negatives observed out of all true negatives expected, and so specificity declines in direct relation to an increase in the number of false positives. For example, a specificity of 0.8 would generate 200 false positives in a sample of 1,000 unlabeled taxa. Specificity, as measured across a wide range in parameters, was predicted to be highest for MW-HR-SIP (1.00 ± 0; ave. ± s.d.) and HR-SIP (1.00 ± 0), substantially lower for qSIP (0.88 ± 0.06), and very low for Heavy-SIP (0.28 ± 0.16) (Figure 4).

Sensitivity is the fraction of true positives observed out of all true positives expected. For example, a sensitivity of 0.7 means that a method failed to detect 30% of the incorporators present. Both qSIP and Heavy-SIP are predicted to have relatively high sensitivity (median values of 0.91 and 0.93, respectively), and the sensitivity of these methods was largely insensitive to the atom % excess of DNA and the number of incorporators (Figure 4). In contrast, the sensitivities of both HR-SIP and MW-HR-SIP were predicted to be highly responsive to the atom % excess of DNA, and the number of incorporators (Figure 4). For these methods, sensitivity is predicted to decline in proportion to the atom % excess ¹³C label in DNA.

Balanced accuracy is calculated as the mean of specificity and sensitivity. The model predicts a tradeoff in balanced accuracy in relation to the atom % excess ¹³C of DNA. MW-HR-SIP had the highest predicted accuracy of any of the 7 methods tested when % atom excess ¹³C exceeded 50%, but qSIP has higher accuracy at lower levels of isotope incorporation (Figure 4). This tradeoff in balanced accuracy resulted from a difference in the tolerance for false positives. For example, MW-HR-SIP produced nearly zero false positives but as a result of its high specificity, it lost sensitivity at lower levels of isotope incorporation. In contrast, qSIP detected labeled taxa across a wider range of isotope incorporation, but it did so at the cost of a large number of false positives.

We also evaluated the effects of sequencing effort and fraction size on DNA-SIP accuracy because SIP experiments often vary in the number of fractions analyzed per gradient and the number of sequences analyzed per fraction. The model shows that sequencing effort has different consequences for method specificity and sensitivity. We found that the specificities of HR-SIP and MW-HR-SIP are predicted to be independent of sequencing effort, while those of Heavy-SIP and qSIP actually declined with sequencing effort (Figure 5). The predicted decrease in specificity for Heavy-SIP and qSIP indicates that the rate of false discovery for these methods increases in proportion to the number of sequences analyzed, while the rate of false discovery in HR-SIP and MW-HR-SIP is unaffected by the number of sequences analyzed. In contrast, we found that method sensitivity improved with sequencing effort for all analytical methods (Figure 5). This predicted increase in sensitivity is caused by an increase in statistical power caused by sampling more sequences from each gradient fraction. We further show that improvements in sensitivity are predicted to be greatest for taxa present at low relative abundance (Figure 6). For example, with MW-HR-SIP the sensitivity of detection for a 50% atom ¹³C enriched OTU present at 0.001 relative abundance is predicted to be nearly zero if less than 1,000 sequences are analyzed per gradient fraction, but sensitivity improves dramatically as sequencing depth increases (Figure 6).

SIP methods also will often vary in the number of fractions collected per gradient, and so we evaluated the effect of fraction size on DNA-SIP accuracy. The model predicts that the use of smaller fractions (i.e., collecting more fractions per gradient) tends to improve sensitivity and overall accuracy for HR-SIP, MW-HR-SIP, and qSIP, though the effects are modest (Figure S7). In contrast, Heavy-SIP is predicted to improve in accuracy when larger fractions are used, though this effect is also somewhat modest (Figure S7).

All DNA-SIP analyses rely upon comparisons made between isotopically enriched experimental treatments and their corresponding unlabeled controls. In real SIP experiments, the composition of replicate post incubation communities are likely to vary somewhat as a result of sample heterogeneity and incubation effects. However, the simulations described above assume random sampling from identical pre-fractionation (post-incubation) community structures. We hypothesized that an increase in variation in community composition between treatment and control samples would decrease the accuracy of DNA-SIP analyses. To test this hypothesis, we generated simulations in which isotope incorporation was held constant (50 atom % excess ¹³C; 10% of OTUs are incorporators) but beta-diversity was varied among 3 replicate treatment and 3 replicate control samples. We varied beta-diversity in two ways: (i) using permutation to vary the rank abundance of a fixed proportion of community members and (ii) varying the proportion of taxa shared between communities. For each simulation scenario, we calculated the mean Bray-Curtis distance among communities in order to provide a real-world metric for gauging the potential accuracy of actual DNA-SIP experiments.

The model predicts that increased beta-diversity among samples impacts the accuracy of DNA-SIP methods (Figure 7). The model predicts that accuracy is impacted more by the number of taxa shared between samples than by differences in taxon abundance (Figure S8). The sensitivity of all methods is predicted to decline as beta-diversity increases, falling from approximately 0.9 for both qSIP and MW-HR-SIP when samples shared 100% of their OTUs to 0.64 and 0.7 for qSIP and MW-HR-SIP, respectively, when samples shared only 80% of their OTUs (Figure S8). Increasing the beta-diversity between samples was predicted to have little effect on the specificity of qSIP but diminished the specificity of HR-SIP and MW-HR-SIP (Figure 7). MW-HR-SIP was predicted to have greater balanced accuracy than qSIP as the Bray-Curtis distance between treatment and control samples increased from of 0.0 to 0.4, but at higher levels of distance the two methods performed with similar accuracy (Figure 7). Regardless, it is clear that sample-to-sample variation has an overall negative impact on DNA-SIP accuracy, and this result emphasizes the importance of minimizing experimental variation between unlabeled controls and labeled treatments in SIP experiments.

So far, we have focused on the accuracy of DNA-SIP methods with respect to the identification of taxa that incorporate isotope into their DNA. However, changes in DNA BD can also be used to quantify the isotope enrichment of DNA from particular taxa. Two approaches have been used to evaluate isotope enrichment from DNA-SIP data: qSIP and ΔBD, with the latter being a complementary analysis to HR-SIP (Pepe-Ranney et al., 2016a). Both ΔBD and qSIP derive quantitative estimates from measuring taxon BD shifts (and thus atom % excess) in the labeled treatment gradient(s) vs. their unlabeled counterparts. The ΔBD method attempts to measure the extent of the BD shift directly from the compositional sequence data, while qSIP utilizes relative abundances transformed by qPCR counts of total SSU rRNA copies. Therefore, ΔBD accuracy likely suffers from compositional effects inherent to HTS datasets, while qSIP accuracy is dependent on qPCR accuracy and variation.

We assessed the quantification accuracy of both methods using the simulations described previously, where either the amount of isotope incorporation or sample beta-diversity was varied. The model predicts that ΔBD produced estimates of isotope incorporation that are closer on average to the true value compared to qSIP, but ΔBD values had much higher variance than qSIP estimates (Figure 8). Furthermore, the predicted variance in ΔBD atom % excess ¹³C estimates increased substantially with even moderate increases in beta-diversity between samples, while the qSIP estimations were largely invariant across the simulation parameter space (Figure S9A). Both methods are predicted to consistently misestimate true ¹³C atom % excess, though the effect was greater for qSIP, with qSIP underestimating ¹³C atom % excess by 30.2–39.2% for fully labeled DNA (Figure 8B).

We further investigated several factors to determine whether they impact the estimation of ¹³C atom % excess from DNA-SIP data. The model predicts that the size of gradient fractions (Figure S10) and the depth of sequencing per fraction (Figure S11) had little impact on estimation of ¹³C atom % excess for OTUs. Furthermore, the model predicts that combining techniques, by using MW-HR-SIP to first identify labeled taxa and then using qSIP to calculate the ¹³C atom % excess of OTUs, reduces the variance of ¹³C atom % excess estimates but does not correct for the systematic underestimation of ¹³C atom % excess (Figures S10, S11). Finally, we attempted to determine why qSIP systematically underestimates the ¹³C atom % excess of OTUs. We find that the degree to which qSIP underestimates ¹³C atom % excess is predicted to increase in proportion to the actual ¹³C atom % excess of each OTU (Figure S12). This outcome could be explained by qSIP's use of weighted averaging in estimating ¹³C atom % excess. The use of weighted averages to estimate the BD of each OTU assumes a roughly Gaussian BD distribution, which we show not to be the case (Figure S1). The BD distributions of “heavy” DNA fragments will be left skewed in a CsCl gradient and the weighted average of a left skewed distribution will cause systematic underestimation of ¹³C atom % excess with the degree of underestimation increasing in proportion to the BD of the DNA, as observed.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.