Stable isotope probing (SIP) facilitates culture-independent identification of active microbial populations within complex ecosystems through isotopic enrichment of nucleic acids. Many DNA-SIP studies rely on 16S rRNA gene sequences to identify active taxa, but connecting these sequences to specific bacterial genomes is often challenging. Here, we describe a standardized laboratory and analysis framework to quantify isotopic enrichment on a per-genome basis using shotgun metagenomics instead of 16S rRNA gene sequencing. To develop this framework, we explored various sample processing and analysis approaches using a designed microbiome where the identity of labeled genomes and their level of isotopic enrichment were experimentally controlled. With this ground truth dataset, we empirically assessed the accuracy of different analytical models for identifying active taxa and examined how sequencing depth impacts the detection of isotopically labeled genomes. We also demonstrate that using synthetic DNA internal standards to measure absolute genome abundances in SIP density fractions improves estimates of isotopic enrichment. In addition, our study illustrates the utility of internal standards to reveal anomalies in sample handling that could negatively impact SIP metagenomic analyses if left undetected. Finally, we present SIPmg, an R package to facilitate the estimation of absolute abundances and perform statistical analyses for identifying labeled genomes within SIP metagenomic data. This experimentally validated analysis framework strengthens the foundation of DNA-SIP metagenomics as a tool for accurately measuring the in situ activity of environmental microbial populations and assessing their genomic potential.


Answering the questions, “who is eating what?” and “who is active?” within complex microbial communities is paramount for our ability to model, predict, and modulate microbiomes for improved human and planetary health. These questions can be pursued using stable isotope probing to track the incorporation of labeled compounds into cellular DNA during microbial growth. However, with traditional stable isotope methods, it is challenging to establish links between an active microorganism’s taxonomic identity and genome composition while providing quantitative estimates of the microorganism’s isotope incorporation rate. Here, we report an experimental and analytical workflow that lays the foundation for improved detection of metabolically active microorganisms and better quantitative estimates of genome-resolved isotope incorporation, which can be used to further refine ecosystem-scale models for carbon and nutrient fluxes within microbiomes.


The explosion of environmental sequencing data in the last decade has fueled a deeper understanding of the role of microbiomes in shaping human health, ecosystem function, and the Earth’s biogeochemical cycles (1). Further advancements in microbiome science require improved experimental approaches that link genomes to their in situ activities. Due to the limitations of culturing techniques, culture-independent methods that reveal in situ functions and link them to taxonomic identities play a crucial role in advancing the field of microbial ecology (2). Stable isotope probing (SIP) is a powerful cultivation-independent tool that links metabolic activity and taxonomic identity of environmental microbes (3, 4). During a DNA-SIP experiment, compounds enriched with heavy stable isotopes (e.g., 13C, 15N, and 18O) are added to the microbial community of interest. The labeled compound is metabolized by active members of the microbial community and incorporated into cellular components, including DNA, during growth (4, 5). As a result, the DNA of these active microbes becomes increasingly isotopically labeled and, therefore, “heavier” compared with the non-labeled DNA from inactive microbes (4, 5). Isotopically labeled DNA, referred to as “labeled” from hereafter, can be physically separated and recovered via isopycnic centrifugation using a CsCl gradient (6). Thus, microbes assimilating labeled compounds in situ can be identified through comparative sequence analysis of the DNA collected at different buoyant densities (BD) along the gradient.
Traditional DNA-SIP studies use 16S rRNA gene sequencing to identify labeled microorganisms (7, 8), and several analysis tools are available for DNA-SIP data (9 - 11). In addition to identifying microbial groups as either labeled or unlabeled, analysis tools such as delta BD (ΔBD) (12) and quantitative SIP (qSIP) (11) can also estimate the extent of isotope assimilation as atom fraction excess (AFE), which is the increase in the isotopic composition of DNA above background levels (11). Measurements of AFE can inform in situ growth rate estimates for specific microbial populations, enabling modeling of microbiome dynamics (13 - 15). Although 16S rRNA gene analyses can taxonomically classify labeled microbes in DNA-SIP studies, the full genomic potential of metabolically active taxa is not always captured due to the difficulty in linking partial 16S rRNA gene sequences to their corresponding genomes (16). Adapting SIP analysis tools for the genomic level, rather than the 16S rRNA gene level, enables genome-centric metagenomic SIP experiments that establish stronger links between genomic information and in situ activity (17).
In recent years, multiple SIP studies have used metagenome sequencing in addition to, or in place of, 16S rRNA gene amplicon sequencing (18 - 23). We refer to this general approach as “SIP metagenomics” from here on to distinguish it from DNA-SIP using 16S rRNA genes. Some recent studies have applied the qSIP approach to shotgun sequencing data to estimate the isotopic enrichment of soil metagenome-assembled genomes (MAGs) (24 - 26). While these represent exciting advancements in the field, SIP metagenomics faces several data analysis and interpretation challenges. For example, estimates of isotopic enrichment depend on accurate measurements of absolute genome abundance, but determining genome abundance from metagenomic data is difficult due to its compositional nature (27 - 30). In addition, outstanding questions remain regarding optimal assembly strategies and the specificity and sensitivity of analysis tools, given varying sequencing depth and genome coverage. Empirically answering these questions requires a defined experiment where the identity of labeled genomes and their level of isotopic enrichment is known a priori. To date, no such empirical study for validating SIP metagenomic sample processing and analysis has been published.
Here, we explored SIP metagenomic sample processing and analysis strategies using an environmental microbiome amended with isotopically labeled Escherichia coli DNA, such that the identity of labeled genomes and their level of enrichment was experimentally controlled. We also investigated the utility of adding internal standards to monitor the quality of density gradient separations and normalize genome coverage levels. With this experimental design, we were able to: (i) compare assembly methods for optimal genome recovery; (ii) determine how sequencing depth and genome coverage influence the detection of labeled genomes; (iii) examine how different approaches for measuring genome abundance impact estimates of AFE; and (iv) compare the sensitivity and specificity of different SIP analysis tools for accurately identifying labeled genomes. Based on our findings, we describe an experimentally validated strategy for SIP metagenomics and provide an R package (SIPmg) that adapts SIP analysis tools for shotgun metagenome sequence data, estimates absolute genome abundance within each fraction using internal standards, and identifies labeled genomes.


To create a ground truth dataset for assessing SIP metagenomics, we generated a microbial community DNA sample where the identity of labeled genomes and their level of enrichment was known a priori (Fig. 1). Specifically, we combined unlabeled DNA extracted from a freshwater pond with aliquots of 13C-labeled E. coli DNA. We created eight levels of E. coli labeling ranging from 0 to 32 atom% 13C enrichment (Table S1 at https://doi.org/10.6084/m9.figshare.22280632). We also added two sets of synthetic DNA oligos at two different stages of sample processing to serve as internal standards (Fig. 1). The six “pre-centrifugation spike-in” standards had different BDs, each reaching maximum abundance in a different and predictable region of the density gradient (Table S2 at https://doi.org/10.6084/m9.figshare.22280632). Deviations from the expected distribution pattern indicated possible problems, such as a disturbance of the density gradient, that might compromise data quality from that sample (Fig. 2). The post-fractionation spike-ins, referred to as “sequins” hereafter (30) (Data Set S1), were added to each fraction after density separation (Fig. 1) to serve as internal calibration standards for calculating absolute genome abundances (Fig. 2). This experimental design provided a controlled dataset for answering questions regarding assembly strategies, genome abundance measurements, the impact of sequencing depth, and the accuracy of various SIP analysis methods.
FIG 1 Experimental design and overview of laboratory steps in the SIP metagenomics workflow. To create a defined SIP experimental sample, DNA extracted from an unlabeled freshwater microbial community was amended with either labeled (13C) or unlabeled (12C) E. coli DNA. Pre-centrifugation spike-ins were added to each sample prior to ultracentrifugation in a CsCl gradient, and post-fractionation spike-ins (sequins) were added to each fraction after density gradient fractionation and collection. These two sets of synthetic DNA oligos served as internal standards to monitor the quality of density separations and normalize genome coverage levels.
FIG 2 The workflow scheme for SIP metagenomic data analysis includes (A) quality filtering of the raw reads and (B) generation of a unique set of medium- and high-quality MAGs used for (C) quantification of absolute taxa abundances and identification of isotope incorporators. The addition of sequins provides the means for calculating absolute bacterial abundances (C, Data Normalization), and pre-centrifugation spike-ins aid in the detection of anomalous samples (C, Outlier Handling).
To develop an empirically validated workflow for SIP metagenomics, we next created the SIPmg R package, which was specifically designed to analyze shotgun sequence data from SIP studies. SIPmg calculates absolute taxon abundances using various methods, such as normalizing relative genome coverage to internal standards (this study) or total DNA concentrations (24, 25). SIPmg provides taxon abundance into the HTS-SIP tool (31) where users can select different methods for identifying isotope incorporators, including qSIP (11), HR-SIP (9), and MW-HR-SIP (10). SIPmg also implements a version of the ΔBD method for estimating isotopic enrichment levels (9). To take advantage of metagenomic data, and similar to Greenlon et al. (25), SIPmg updates the qSIP model to use the observed GC content of assembled genomes rather than the estimated GC content used in qSIP analysis of 16S rRNA gene data (11). Finally, to correct for multiple comparisons, i.e., testing for significant isotope enrichment in multiple MAGs, SIPmg can adjust the confidence intervals around bootstrapped estimates of AFE using a variation of false discovery rate correction (32). With the SIPmg package, we evaluated the performance of different analysis approaches using our ground truth SIP metagenomics dataset.

Maximizing recovery of metagenome-assembled genomes using individual and combined assemblies

In contrast to a typical metagenome sample, community DNA in an SIP experiment is separated into multiple fractions based on BD prior to sequencing (Fig. 1). Differences in GC content and levels of isotopic enrichment result in a non-random distribution of microbial genomes across the density gradient, and sequencing each density fraction provides multiple options for assembly and binning. To determine the optimal strategy for maximizing MAG recovery, we compared assembly of the intact unfractionated sample, separate assemblies of each individual fraction, co-assemblies of all fractions derived from the same initial samples, and a massive combined assembly using MetaHipMer (33) of all fractions from all samples. The latter three strategies all used the same 1,418 Gbp of sequence data from hundreds of sequencing libraries and grouped in different ways for each strategy, while the unfractionated assembly used only 47 Gbp from one sequencing library. Each assembly was then independently binned using MetaBAT2 (34).
A total of 2,022 MAGs were generated across all assemblies, of which 248 were high quality, 447 were medium quality, and 1,327 were low quality as defined by the minimum information about metagenome-assembled genomes (MIMAG) reporting standards (35) (Data Set S2). The MetaHipMer assembly produced more MAGs than any other strategy. A total of 235 MAGs were recovered from the MetaHipMer assembly, of which 136 were medium- or high quality (Fig. 3A). Estimates of average MAG completeness and contamination for each assembly type were not substantially different (Fig. S1).
FIG 3 Comparison of metagenome assembly approaches for the SIP metagenome dataset generated from spiking E. coli into unlabeled DNA from a freshwater microbiome. (A) Average number of medium- and high-quality MAGs recovered from different assembly approaches. (B) Venn diagram showing the number of unique and shared MAG clusters. (C) Compositional differences at the class level recovered from different types of assemblies (I—intact metagenome assembly with MetaSPAdes, F—separate fractions assembled with metaSPAdes [n = 371 assemblies], S—all fractions within each replicate co-assembled with metaSPAdes [co-assembly of all fractions sequenced for a single SIP replicate sample, n = 24 assemblies], M—combined assembly of all fractions using MetaHipMer; for F and S the average number of MAGs was calculated, whiskers represent standard deviation across assembly type).
Next, we deduplicated all the medium- and high-quality MAGs recovered from all assemblies to determine whether any approach generated unique MAGs that were not present in other assembly types (Fig. 2B). We first grouped MAGs with average nucleotide identities of ≥96.5 and alignment fractions of ≥30% into a total of 148 unique clusters (36), then selected a single representative MAG for each cluster. Of these, 120 MAG clusters were exclusively produced by MetaHipMer. Twelve MAG clusters did not include any MetaHipMer-generated MAGs, and 11 of these clusters contained at least one MAG generated from the assemblies of individual fractions (Fig. 3B). Assembly of the intact unfractionated microbiome did not produce any unique MAGs (Fig. 3B), presumably because the sequencing depth for the unfractionated sample (47 Gbp) was much smaller than the total sequencing depth of all the individual fraction assemblies, the sample-wise combined assemblies, and the MetaHipMer assembly (1,418 Gbp). The different assembly strategies also produced MAGs with different taxonomic compositions. For example, MAGs derived from the MetaHipMer assembly accounted for an additional nine classes that were not present in other assemblies (e.g., Anaerolineae, Andersenbacteria, Babeliae, Chlamydia, among others) (Fig. 3C). Most MAGs that were unique to the MetaHipMer co-assembly had lower coverage than MAGs recovered by other assembly approaches (Fig. S1). This suggests the MetaHipMer co-assembly captured more of the lower abundance MAGs in the samples than other assembly approaches, possibly due to the higher coverage levels that resulted from combining reads from all libraries into one assembly (33). These results indicate that employing multiple assembly strategies and dereplicating the resulting MAGs can maximize genome recovery in SIP metagenomics studies.

Anomalous sample detection using pre-centrifugation spike-in controls

As part of the quality control process, we devised an approach for detecting anomalous samples whose pre-centrifugation spike-in sequences displayed aberrant distributions along the BD gradient (Fig. 2C). We added six synthetic spike-ins to our samples prior to ultracentrifugation, and each spike-in had a different density based either on its GC content or the artificial introduction of 13C-labeled nucleotides during oligo synthesis (Table S2 at https://doi.org/10.6084/m9.figshare.22280632); therefore, each spike-in has a distinct and predictable peak in coverage along the BD gradient. Deviations from the expected spike-in distribution patterns may indicate events such as cross-contamination, library misidentification, or accidental disturbances of the density gradient significant enough to distort the distribution of MAGs throughout the gradient, all of which would introduce error into the downstream analysis. We identified three biological replicates with anomalous spike-in distribution patterns (Fig. S2), and these samples were removed from downstream analyses to avoid the introduction of extraneous noise. This example illustrates the utility of internal standards to illuminate quality control problems in SIP experiments that would otherwise go undetected.

Normalizing genome coverage to quantify DNA isotope incorporation

Accurate abundance measurements are critical for determining levels of isotopic labeling. Briefly, models such as qSIP and ΔBD estimate a taxon’s AFE based on differences between its weighted BD in unlabeled controls and isotope-amended treatments (9, 11, 37), where weighted BD is calculated from the taxon’s abundance within each density fraction (see Materials and Methods, Equations 5 and 6). For amplicon-based qSIP studies, the relative abundance of a taxon is normalized to the total number of 16S rRNA gene sequences within each fraction determined by qPCR (11). Estimating abundance in SIP metagenomic studies is more complicated, since shotgun sequencing lacks an equivalent method to 16S rRNA gene qPCR for absolute abundance scaling. Previous SIP metagenomic studies multiplied relative genome coverage with the total DNA concentration of each fraction (25, 26), which is a reasonable approach, although it does not account for potential variability introduced during DNA recovery, library creation, and sequencing of each fraction (29, 30, 38). By adding sequins to density fractions before DNA precipitation and recovery, we explored an alternative normalization strategy for measuring absolute abundance that could also account for variability in the downstream processing steps (24). In this approach, genome coverage within each fraction can be converted into absolute abundances through normalization based on the known concentration and observed coverage of the sequin internal standards. The AFE of each genome can then be estimated from these abundance measurements.
Our experimental design, where isotopic enrichment levels were known a priori, provided an opportunity to compare different approaches for calculating genome abundances and determine their impact on estimates of taxon AFE values (Table 1; Fig. S3). More specifically, we compared the expected AFE values for labeled E. coli with AFE estimates from the qSIP model made with different approaches for calculating abundance, including absolute abundance derived from normalization to sequins (Fig. 4A), absolute abundance estimated by multiplying either relative abundance or relative coverage with total DNA concentration (Fig. 4B and C, respectively), and relative coverage without conversion to absolute abundance (Fig. 4D). Results from all of the abundance normalization strategies we tested are provided in Fig. S3 and Table S3 (https://doi.org/10.6084/m9.figshare.22280632). Any genome other than E. coli that was identified as labeled was considered a false positive, whereas failure to identify E. coli as labeled was considered a false negative.
TABLE 1 Performance of different approaches for calculating genome abundance across density fractions based on the results from spiking 13C-labeled E. coli DNA into background DNA of an unlabeled freshwater communitya
MethodProcedureSpecificitySensitivitySpearman correlation between estimated and true AFE (P-values)
Absolute abundance using sequinsRegression using sequin coverage and concentration0.9930.8570.86 (0.012)
Absolute abundance using total DNA concentrationProduct of relative abundance and DNA concentration (25)0.9910.7140.83 (0.021)
Product of relative coverage and DNA concentration (24)0.9220.5710.38 (0.4)
Relative coverageRelative coverage of MAGs in each fraction0.9990.5710.77 (0.041)
AFE was predicted using the qSIP model. Specificity was estimated as (true negatives)/(false positives + true negatives). Sensitivity was estimated as (true positives)/(true positives + false negatives).
FIG 4 Comparison of predicted AFE versus the expected AFE of E. coli using different approaches for measuring genome abundance across the density gradient. The qSIP method was used to estimate AFE in all cases. Genome abundance in each density fraction was determined by (A) normalization to sequin internal standards, (B) multiplying relative abundance with DNA concentration following Greenlon et al. (25), (C) multiplying relative coverage with DNA concentration following Starr et al. (24), and (D) relative coverage without additional normalization. For all comparisons, please refer to Table S3 (https://doi.org/10.6084/m9.figshare.22280632). Error bars represent the standard deviation of AFE calculated using the qSIP method’s bootstrapping approach. The expected AFE for each condition is in parentheses, and additional details about conditions, including replicate numbers, are provided in Table S1 (https://doi.org/10.6084/m9.figshare.22280632). pcor and preg correspond to the P-values for the Spearman correlation and the linear regression F-statistic, respectively. The intercepts determined by linear regression were not significantly different from zero (P-value > 0.05) in any method for estimating abundance.
Abundance estimates derived from the sequin approach outperformed all other approaches based on combinatorial assessment of specificity (lower false positives), sensitivity (lower false negatives), and the Spearman correlation between expected and predicted AFE values (Fig. 4 and Table 1; Table S3 at https://doi.org/10.6084/m9.figshare.22280632). The two approaches using total DNA concentrations underestimated levels of AFE, and one approach did not produce statistically significant linear regressions (P-value >0.05) between expected and estimated AFEs (Fig. 4B and C; Table S3 at https://doi.org/10.6084/m9.figshare.22280632), although the sensitivity for detecting labeled E. coli was the same or better than sensitivity using relative coverage (Table 1). Relative coverage produced the highest specificity, although it had lower sensitivity than the normalization approach using sequins (Fig. 4D and Table S3 at https://doi.org/10.6084/m9.figshare.22280632). These results suggest that internal quantification standards can improve estimates of genome abundance and AFE.

Comparison of various SIP analysis methods

In addition to qSIP, other analysis methods such as ΔBD (9), HR-SIP (9), and MW-HR-SIP (10) can identify labeled taxa. We compared all four approaches for their ability to accurately identify isotope incorporators in our defined SIP metagenomic dataset. We also compared estimates of E. coli AFE predicted with the ΔBD and qSIP methods; HR-SIP and MW-HR-SIP do not provide quantitative estimates of enrichment. For all methods, absolute genome abundances were determined by normalization to sequins.
Estimates of E. coli AFE made with the qSIP model more closely matched expected isotopic enrichment levels than did estimates from the ΔBD method (Fig. 5). The qSIP approach also had higher specificity than the ΔBD method, producing only 7 false positives across all conditions compared with 12 false positives, respectively (Table 2). The MW-HR-SIP approach had the fewest false positives of any method, with only four across all conditions, while maintaining the same sensitivity, in terms of false negatives, as the qSIP method (Table 2). The sensitivity and specificity of HR-SIP were lower than both MW-HR-SIP and qSIP methods (Table 2). Based on these results, we selected qSIP and MW-HR-SIP for further evaluation.
FIG 5 Comparison of AFE estimates produced by the (A) qSIP and (B) ΔBD methods using the amended metagenome where levels of E. coli isotopic enrichment were known a priori. Both of these methods used sequin-based normalization for estimating genome abundance. Error bars represent the standard deviation of AFE calculated using the qSIP method’s bootstrapping approach. The expected AFE of E. coli within each treatment condition is given in parentheses. preg and pcor correspond to the P-values for the linear regression and Spearman correlation, respectively. The intercepts determined by linear regression for qSIP and ΔBD models were not significantly different from zero (P-value > 0.05).
TABLE 2 Comparison of methods to identify isotopically labeled genomesa
Incorporator identification methodFalse positivesSpecificitySensitivityBalanced accuracy
qSIP model70.9930.8570.925
ΔBD method120.9840.8570.921
Evaluations were based on absolute genome abundances obtained by normalizing coverage to internal sequin standards using the sequin approach. Specificity and sensitivity were averaged across the seven treatment conditions.

Lower limits of genome coverage for reliable detection of isotope labeling

Next, we examined how sequencing depth affected our ability to detect isotope incorporation. As demonstrated above, the accuracy of abundance measurements impacts the accuracy of AFE estimates, and these abundance measurements are derived from genome sequence coverage. The relative abundance of microbial taxa comprising complex communities can vary by orders of magnitude; thus, genome coverage within sequencing libraries can vary similarly (39). This suggests that AFE estimates might be less reliable for taxa with low coverage. To determine the lowest depth of coverage at which an AFE could be accurately estimated, we performed qSIP and MW-HR-SIP analyses after subsampling E. coli reads to 10%, 1%, 0.1%, 0.01%, and 0.001% of their initial levels (Table S4 at https://doi.org/10.6084/m9.figshare.22280632). In the respective subsampled datasets, E. coli had an average total coverage ranging from 0.01× to 1,400× coverage. Here, “total coverage” refers to the cumulative coverage across all density fractions of an individual sample.
The qSIP model consistently identified E. coli as labeled when mean total coverage was ≥1× (Table S5 at https://doi.org/10.6084/m9.figshare.22280632). The correlation coefficient between actual and predicted AFEs was 0.83 within this coverage range (P-value < 0.05; Fig. S4 and Table S6 at https://doi.org/10.6084/m9.figshare.22280632). However, at total coverages <1×, qSIP failed to detect E. coli as labeled in several experimental conditions (Fig. S4 and Table S6 at https://doi.org/10.6084/m9.figshare.22280632). The MW-HR-SIP method was also less sensitive at lower coverage levels, and at 100× mean total coverage, it only detected E. coli as labeled in three out of seven experimental conditions (Table S5 at https://doi.org/10.6084/m9.figshare.22280632). These data suggest that estimates of isotope enrichment are less reliable, in general, when genome coverage is low.

Strategies to improve accuracy of detecting isotopically labeled genomes

To improve the accuracy of SIP metagenomic analysis, we explored different strategies to reduce the number of genomes incorrectly identified as labeled (i.e., false positives). For example, the number of false negatives increased as coverage decreased; therefore, we tested whether implementing minimum genome coverage requirements could reduce the number of false positives. Excluding genomes with mean total coverages < 10× reduced the total number of MAGs analyzed from 147 to 113 and reduced false positives identified by qSIP from 7 to 4 without increasing false negatives (Table S7). This improved the balanced accuracy from 0.925 to 0.927. Raising the minimum mean total coverage to 17× eliminated all false positives (Fig. S5), but this also reduced the number of remaining MAGs analyzed to only 68. We did not test coverage limits for MW-HR-SIP because the method struggled to detect E. coli as labeled when coverage dropped below 100× (Table S5 at https://doi.org/10.6084/m9.figshare.22280632) and applying a threshold of 100× would have limited our analysis to only 17 genomes (Table S7 at https://doi.org/10.6084/m9.figshare.22280632). These results suggest that excluding genomes with low coverage can decrease false positives and increase balanced accuracy. Although the definition of “low coverage” will vary based on experimental conditions and individual assessments of the trade-offs between sensitivity and specificity, these results also suggest that confidence in the identification of labeled genomes should decrease along with their coverage levels.
We also investigated if false positives could be reduced by implementing a minimum level of isotopic enrichment necessary for a genome to be considered labeled. That is, rather than simply requiring genomes to be significantly greater than 0% AFE, which is the default setting of the qSIP approach (11), we examined different minimum AFE thresholds ranging from 2% to 12.5% (Table S8 at https://doi.org/10.6084/m9.figshare.22280632). A genome was considered to be labeled if the lower bound of its AFE 95% CI was greater than the minimum AFE threshold. With AFE thresholds between 2% and 6%, total false positives dropped from 7 to 3 across all experimental treatments, but E. coli was no longer identified as labeled in one experimental condition. The balanced accuracy was also reduced from 0.925 without AFE thresholds to 0.856 with a 6% AFE threshold (Table S8 at https://doi.org/10.6084/m9.figshare.22280632). False positives were completely eliminated with a minimum AFE threshold of 12.5%, but sensitivity was so poor (0.286) that E. coli was only identified as labeled in 2 out of 7 conditions (Table S8 at https://doi.org/10.6084/m9.figshare.22280632). Minimum AFE limits could not be tested with MW-HR-SIP analysis because this method does not estimate levels of isotopic enrichment. Together, these results illustrate a trade-off between sensitivity and specificity when increasing the minimum AFE threshold above zero and suggest that false positives can be reduced by increasing the AFE threshold at the potential cost of losing sensitivity for the detection of minimally labeled taxa.
The number and identity of false positives varied across SIP analysis methods, presumably due to differences in the methods’ underlying algorithms. Therefore, we hypothesized that the number of false positives might be reduced by taking the consensus of different analysis methods, i.e., requiring that two separate models predict a MAG is labeled. All false-positive MAGs found in qSIP analysis were also false positives in ΔBD analysis, thus taking the consensus of these two methods did not produce fewer false positives than qSIP alone (Table S9 at https://doi.org/10.6084/m9.figshare.22280632). In contrast, there was no overlap in the identity of false-positive MAGs between the qSIP and MW-HR-SIP methods, and a union of their results completely eliminated false positives without producing any false negatives (Table S9 at https://doi.org/10.6084/m9.figshare.22280632). However, we found it more advantageous to apply MW-HR-SIP and qSIP sequentially rather than independently. MW-HR-SIP had greater specificity than qSIP; therefore, it was used as a first-pass filter to detect putatively labeled genomes while minimizing false positives. This subset of putatively labeled genomes was then re-analyzed with the qSIP model. Only genomes first identified as labeled by MW-HR-SIP and later confirmed with a significantly positive AFE by qSIP were labeled. Applying the tools in series reduced the number of multiple hypotheses tested (i.e., MAGs tested for enrichment), which subsequently increased the statistical power for AFE estimation. That is, without the initial reduction in identified incorporators, the qSIP analysis would have otherwise included all MAGs in its statistical comparisons between treatment groups, resulting in a smaller P-value required for significance with multiple hypothesis testing. The increased statistical power obtained by applying the models in series resulted in tighter confidence intervals for the AFEs of E. coli (Table S10 at https://doi.org/10.6084/m9.figshare.22280632). These results indicate that using a combination of analysis tools can reduce false-positive detection, although the tools used and their order of application may vary depending on preferences for sensitivity versus specificity.


DNA-SIP has been an established method in microbial ecology for many years and has primarily relied on 16S rRNA gene sequencing to identify active taxa (4, 5, 11, 15, 18, 40). With decreases in sequencing costs and increases in compute capacity, DNA-SIP studies can now utilize shotgun metagenomic sequencing to establish links between population genomes and in situ activities (24 - 41 - 43 - 43). In addition, automated sample preparation substantially increases the potential scale of SIP metagenomic studies and allows for more biological replication (26). However, the growth of SIP metagenomics also depends on adapting analysis tools to work with shotgun metagenomic data and validating their performance. To this end, we investigated an artificial SIP metagenome that enabled empirical testing of sample processing and data analysis strategies. Our results suggest some potential improvements to SIP metagenomic methodologies that can serve as a foundation for future advances.
Comparing assembly strategies for SIP metagenomic data was a key goal of our study. Previous SIP studies have used different strategies, including assembling unfractionated DNA, assembling individual SIP fractions, and co-assembling several fractions (24 - 44 - 45 - 45). However, it was not clear which assembly strategy produces the most medium- and high-quality MAGs. For instance, in computationally simulated SIP experiments, the co-assembly of multiple fractions improved MAG recovery compared with the assembly of unfractionated DNA (45). In addition, the large amount of sequence data used in co-assemblies can recover rare genomes that would otherwise be lost due to insufficient coverage in smaller assemblies of individual datasets (33). Conversely, individual assemblies can outperform co-assemblies in samples where high levels of microdiversity impede contig formation (46 - 48). Here, we found that co-assembly of all density fractions generated the most medium- and high-quality MAGs, which agrees with two recent SIP metagenomics studies (25, 26). However, we also found that merging binning results from individual fraction assemblies and larger co-assemblies via MAG de-replication provided more medium- and high-quality MAGs than did co-assembly alone. We posit that this approach reaps the benefits of both strategies: it provides higher read recruitment for assembling rare genomes in co-assemblies and also leverages lower microdiversity in individual fraction assemblies. Optimal assembly strategies may differ for other environmental samples, and these strategies must be re-evaluated as sequencing and assembly methods evolve, but our results suggest that SIP metagenomic studies can benefit from employing multiple assembly approaches to maximize genome recovery.
Processing DNA-SIP samples is laborious, but semi-automated protocols simplify lab work and enable high-throughput SIP metagenomic studies (26). Indeed, increasing the number of biological replicates and sequencing more density fractions per replicate can improve the detection of labeled taxa (41). However, the opportunities for accidental mistakes, such as cross-contamination, sample mix-ups, or clerical errors, also increase when processing dozens of samples and hundreds of density fractions. In addition, slight mishandling of ultracentrifuge tubes can disturb delicate CsCl gradients (8) and potentially alter genome distributions along the gradient. If left undetected, these types of accidents could produce inaccurate weighted BD estimates, adding extra noise to the data analysis and even compromising results. In this study, we found that including pre-centrifugation spike-ins, which had distinct and predictable distribution patterns along the gradient, helped us identify and remove problematic samples before they negatively impacted our analyses. Including internal standards can mitigate potential errors and enhance the quality of large complex SIP studies with many replicates. Moreover, with careful design and additional development, internal standards might someday correct for variability introduced during sample processing (41) instead of simply flagging samples for removal. Internal standards can be easily incorporated into automated SIP metagenomics protocols (26), where they can improve the quality of SIP metagenomic results, and if adopted broadly, potentially serve as consistent fiducial reference points that facilitate inter-comparisons of different SIP studies.
Accurate measurements of genome abundance along the BD gradient are essential for identifying labeled genomes and determining their level of isotopic enrichment (11). However, the compositional nature of metagenomic data, and the variability introduced during sample processing and sequencing, can hamper quantitative estimates of genome abundance (27 - 30 - 49). Internal quantification standards can mitigate process variability and provide absolute abundance estimates of genes, transcripts, and genomes from metagenome and metatranscriptome data (30, 38, 50 - 53). Based on these findings, we hypothesized that adding internal standards to density fractions (sequins) could improve abundance measurements, thereby improving isotope enrichment measurements. Indeed, estimates of AFE in our study were more accurate using absolute abundances derived from sequin normalization compared with AFE estimates using other strategies.
Multiple factors could explain the more accurate estimates of isotopic labeling enabled by internal quantification standards. For one, sequins may have mitigated any variation introduced during library creation and sequencing (30). Additionally, sequins may have corrected for differences in DNA recovery among fractions that would have otherwise gone unnoticed and negatively impacted abundance measurements. That is, after collecting CsCl fractions, each fraction separately undergoes PEG precipitation and desalting before DNA concentrations are measured (26). Absolute abundances calculated using DNA concentrations assume identical DNA recovery efficiencies (24, 25), so any stochastic or systematic variability in the percent of DNA recovered would lead to errors in absolute abundance measurements. Conversely, sequins track and mitigate variability in DNA recovery when they are added to fractions before the desalting steps, as was performed here. Therefore, if DNA recovery efficiency varied among fractions, then we would expect absolute abundances derived from sequins to be more accurate than estimates derived from DNA concentration measurements. Without internal standards, variability introduced during DNA recovery, library construction, and sequencing is unknowingly propagated as noise into downstream SIP analyses. This undetected variability can potentially lead to errors that impact predictions of isotope enrichment.
The various SIP analysis methods examined in this study use different approaches to detect labeled microorganisms, and these differences could impact the sensitivity and specificity of their predictions. The accuracy of different SIP analysis methods has not to our knowledge been assessed with metagenomic data until now, but in silico simulations of 16S rRNA gene-based SIP data revealed that MW-HR-SIP had higher balanced accuracy than the other analysis methods (31). The qSIP model also generated more accurate AFE estimates than the ΔBD method in those simulations. We observed similar patterns by comparing analysis methods using our experimentally designed SIP microbiome. In addition, we found that the consensus of multiple approaches, i.e., MW-HR-SIP and qSIP, produced higher accuracy results than any single method alone. Future SIP metagenomic studies might increase confidence in identifying isotope-incorporating taxa by employing these two independent strategies, although the higher confidence in true positives might come at the cost of missing labeled genomes with lower coverage. Regardless of the analysis tools used, analyzing more biological replicates is another simple strategy to increase accuracy (41). As SIP analysis methods evolve, reassessing their performance with deeper sequencing, more replicates, and an improved artificial SIP microbiome (e.g., more species at different AFE levels) will provide additional insights into their accuracy and limitations.
Altogether, we used an environmental metagenome amended with isotopically labeled E. coli DNA to assess the performance of different analysis approaches and established an experimentally validated workflow for SIP metagenomics. The “wet-lab” aspects of the workflow include the addition of pre-centrifugation spike-ins for quality control and post-fractionation sequins for genome quantitation along the BD gradient. The “dry-lab” aspects entail an absolute genome normalization in each density fraction, and a modified qSIP model tailored to handle genome-resolved metagenomic datasets to calculate AFE. We also explored strategies to more accurately identify isotope incorporators, such as limiting analysis to taxa with coverage and isotope enrichment levels above minimal thresholds and using the consensus of multiple SIP analysis tools to detect labeling using our newly developed SIPmg package. These additional strategies hold promise for improving the accuracy of SIP metagenomic results, although the specifics of how and when to apply them will depend on the study design and individual preferences regarding the trade-offs between specificity and sensitivity. We believe this validated analysis workflow will increase the reliability of SIP metagenomic findings, enable standardization across studies, and facilitate the use of SIP data in modeling microbially mediated processes.


DNA collection and microbiome amendments

To create a microbiome where the identity of labeled genomes and their level of enrichment was known a priori, we first extracted DNA from a bacterial isolate grown in 13C-labeled glucose. E. coli K-12 wild-type cells were grown in M9 minimal salts media (Teknova; M8005). Glucose was added at a final concentration of 20 mM and was the sole carbon source. DNA with different levels of 13C enrichment was produced by varying the ratio of unlabeled glucose to uniformly labeled 13C₆-D-glucose (Cambridge Isotope Laboratories; CLM-1396; 99 atom %), e.g., DNA extracted from cultures grown in a ratio of 4:1 of unlabeled:labeled glucose was expected to have an enrichment of approximately 20 atom%. Cultures grown overnight in LB were transferred into labeled media at 5,000-fold dilution (i.e., 2 µL into 10 mL labeled media), grown at 37°C, and harvested at midlog phase. DNA was extracted using the Wizard genomic DNA purification kit (Promega; A1120) and quantified using the QuantIT dsDNA high sensitivity assay kit (ThermoFisher; Q33120). To verify the level of 13C enrichment, E. coli DNA was hydrolyzed into nucleosides before analysis by liquid chromatography-mass spectrometry (LC-MS). Briefly, using a modification of the approach described previously (54), 25 ng of DNA was hydrolyzed by incubating for 2 h at 37°C in a 61-µL cocktail containing a 1× final concentration of DNase I reaction buffer (NEB; M0303S) and four units each of DNase I (NEB; M0303S), alkaline phosphatase (Takara; 2120A), and phosphodiesterase (MPbio; 100978). The intensity of each deoxyadenosine isotopologue, ranging from 0 to 10 13C atoms, was measured by LC-MS to calculate AFE (55). Isotopic enrichment levels of E. coli DNA are reported in Table S1 (https://doi.org/10.6084/m9.figshare.22280632).
DNA from a complex microbial community was recovered from an outdoor, man-made pond located at the Joint Genome Institute. Pond water was pre-filtered through a 5-µm mesh before collection onto 0.2-µm Supor filters (Pall; 47mm diameter). DNA was extracted from filters using a DNeasy PowerWater kit (Qiagen; 14900-50-NF).
Replicate samples were prepared for ultracentrifugation by combining 900 ng of microbiome DNA with 50 ng of E. coli DNA. For samples with isotopically labeled DNA, the ratio of unlabeled-to-labeled DNA was adjusted to produce an average overall level of enrichment, e.g., 20 ng of unlabeled E. coli DNA combined with 30 ng of 54% AFE E. coli DNA would have an overall AFE of 32%. The specific ratios of unlabeled:labeled DNA for different conditions, and the number of replicates per condition, are described in Table S1 (https://doi.org/10.6084/m9.figshare.22280632). These mixtures of unlabeled and labeled E. coli DNA produced distribution patterns along the density gradient that differ from those that would have been produced from uniformly labeled DNA with the same overall AFE. However, the SIP analysis tools tested in this study estimate AFE based on differences in weighted BD measurements between controls and treatments (see Equations 1 and 6 below), and the weighted BDs for mixtures are the same as those for uniformly labeled DNA with the same overall AFE (see Equation 5 below). Thus, for the purpose of evaluating SIP analysis tools, mixtures of unlabeled and labeled DNA can be used to simulate various levels of AFE.

Synthetic pre-centrifugation DNA spike-ins

A set of six synthetic DNA fragments were added to mixtures of DNA from isolates and the complex microbiome to track the ultracentrifugation and fraction collection steps. These fragments were approximately 2 kbp in length with GC content of 37%–63% (Table S2 at https://doi.org/10.6084/m9.figshare.22280632). To change the distribution of fragments across the density gradient, some fragments were artificially enriched with 13C through PCR by adjusting the ratio of unlabeled dNTPs and uniformly labeled 13C dNTPs (Silantes Gmhb; 120106100; >98 atom%) (Table S2 at https://doi.org/10.6084/m9.figshare.22280632). Briefly, DNA was amplified for 30 cycles by adding 0.5-µL Phusion high-fidelity (HF) DNA polymerase (NEB; M0530S), 10 µL of 5× Phusion HF Buffer, 1 µL of 10 mM dNTPs (final conc. labeled/unlabeled mixture), 2.5 µL each 10 µM Forward and Reverse Primer, and 31.5 µL of nuclease-free water. PCR products were purified using AMPure XP beads (Beckman Coulter; 63880) and pooled in equimolar ratios to create a set of pre-centrifugation DNA spike-ins. These pre-centrifugation spike-ins were added at 1% by mass of the DNA mixture, e.g., 10 ng of synthetic fragment pool added to 1 µg of microbial DNA mixture.

Gradient separation, sequin addition, and fraction purification

Following Nuccio and colleagues (26), samples were centrifuged at 44,000 RPM (190,600 g) for 120 h at 20°C in a VTi 65.2 Rotor (Beckman Coulter; 362754). For each sample, 24 fractions of 220 µL were collected into a 96-well plate using an Agilent 1260 fraction collector running at flow rate 250 µL/min while using mineral oil as the displacement fluid. Fraction density was determined using a Reichert AR200 refractometer.
Before purifying DNA from CsCl fractions, an additional set of 80 synthetic DNA fragments, or sequins (30), were added to each fraction as an internal standard for subsequent quantitative metagenomic analysis. Lyophilized pellets of sequins were obtained from the Garvan Institute of Medical Research (https://www.sequinstandards.com). Pellets were resuspended in TE Buffer (10 mM Tris, 0.1 mM EDTA, pH 8.0), and the concentration was measured with QuantIT dsDNA high-sensitivity assay kit (ThermoFisher; Q33120). Of the 24 BD fractions collected for each sample, we selected 16 to move forward with library creation and sequencing based on the range of BD they spanned. These 16 fractions were amended with sequins. To compensate for expected differences in the amount of DNA recovered from different densities, the middle eight fractions received 300 pg of sequins while the four fractions on either tail received 100 pg of sequins.
After sequin addition, DNA was recovered by adding a 250-µL solution of 36% polyethylene glycol (PEG) 6000 and 1.6 M NaCl to each fraction and incubating overnight in 4°C. Plates were centrifuged at 3,214 × g for 1.5 h at 20°C to pellet DNA. Pellets were washed with 300 µL of 70% chilled ethanol, centrifuged at 3,214 × g for 45 min at 20°C, and resuspended in 30 µL of TE Buffer (10 mM Tris, 0.1 mM EDTA, pH 8.0). Purified DNA was quantified using Quant-IT dsDNA high-sensitivity assay kit (ThermoFisher; Q33120).
Sequins were added to each fraction before PEG precipitation and DNA quantification steps; therefore, the amount added was based on the expected sample DNA concentrations. Tailoring sequin additions to actual sample DNA concentrations, as opposed to estimates, is preferable to ensure optimal coverage in sequencing data. After completing analysis of the amended microbiome, we sought to improve sequin additions for future studies by measuring DNA levels before PEG precipitation when DNA was still in concentrated CsCl. Additional details are provided in the Supplementary Information.

Library creation and sequencing

Sequencing libraries were generated from the 16 middle fractions of each sample using Nextera XT v2 chemistry (Illumina) with 12 PCR cycles. Concentrations and size distributions of each library were determined on a Fragment Analyzer (Agilent). Libraries were pooled at equal molar concentrations within the range of 400–800 bp, and the pool was size selected to 400–800 bp using a Pippin Prep 1.5% agarose, dye-free, internal marker gel cassette (Sage Science). For each library, 2 × 150 bp paired-end sequencing was performed on the Illumina Novaseq platform using S4 flowcells (Table S6 at https://doi.org/10.6084/m9.figshare.22280632).

Metagenome assembly and binning

Raw reads were filtered and trimmed using RQCFilter2 software according to the standard JGI procedures (https://jgi.doe.gov/data-and-tools/software-tools/bbtools/bb-tools-user-guide/data-preprocessing/). Then, one of the four strategies was used to perform contigs assemblies: (i) an assembly of unfractionated SIP sample with metaSPAdes(v3.15.2) (56); (ii) a single fraction assembly with metaSPAdes (371 assemblies); (iii) a single sample co-assembly with metaSPAdes (co-assembly of all fractions sequenced for a single SIP replicate sample, 24 assemblies); and (iv) an experiment-wise co-assembly with MetaHipMer(v2.0.1.2) (assembly of all fractions across all replicates) (33). Assembly and genome mapping parameters are reported in the Supplementary Methods. We generated 397 assemblies in total. Quality assessment metrics for each assembly were calculated using QUAST(v5.0.2) (MetaQUAST mode) (Data Set S3 ) (57). Each assembly was then independently binned with MetaBAT(v2.12.1) (58). For each generated MAG, we used GTDB-Tk(v2.0.0) (GTDB R95) (59) to assign a taxonomic classification. To assess the quality of MAGs, we used CheckM(v1.1.3) (60) and QUAST(v5.0.2) (61). The MetaHipMer combined assembly was annotated using the JGI metagenome annotation workflow (58) and is available through IMG/M (62) under taxon identifier 3300045762.

MAG deduplication and mean scaffold coverage calculations

Medium- and high-quality MAGs recovered from all assembly strategies were deduplicated to remove redundant versions of each draft genome (35). The genome-wide average nucleotide identity (gANI) and the alignment fraction (AF) were calculated for each possible MAG pairwise comparison (36). Next, the lowest pairwise values of gANI and AF were utilized for each MAG comparison, followed by clustering using single linkage to group MAGs based on species-level delineations (e.g., gANI ≥96.5 and AF ≥30) as defined by Varghese and colleagues (36). MAGs that did not cluster with other MAGs were considered singletons. Following clustering, we used completeness, contamination, and total length values to select a single representative MAG for each cluster. Sequences of all spike-ins and sequins were concatenated with the final set of MAG contigs, and this contig set was then used as a reference for read mapping across all density fractions (see Supplementary Methods). The average contig coverage of MAGs, spike-ins, and sequins in each fraction was calculated and used in the downstream analysis.

Quality control of SIP data using pre-centrifugation spike-ins

Before performing SIP analysis, we first removed mishandled samples from our dataset. For this purpose, we identified the peak of absolute concentration distributions across the density gradient for each labeled pre-centrifugation spike-in. If the spike-in distribution patterns did not match the expected order along the density based on the theoretical estimated density of the spike-in (given its GC content and C13/C12 ratio), then the sample was considered potentially problematic and removed from the analysis.

Estimating the absolute abundance of MAGs across density fractions

To determine the extent of isotope incorporation into genomes, it is first necessary to measure genome abundance across the density gradient. We explored several ways to measure genome abundance in the SIP dataset, which are implemented as part of the SIPmg R package (see Code Availability).
First, we obtained absolute concentrations of genomes across the density gradient using the approach proposed by Hardwick and colleagues (30), in which sequins were used as internal reference standards to scale coverages into absolute concentrations. Briefly, the average MAG coverage within a given fraction (metagenome) was scaled into units of molarity using regression analysis based on known molarity of 80 sequins and their average coverages. Molar concentrations of the sequins in the added standard mixture were obtained from the manufacturer (Garvan Institute of Medical Research). For regression analyses, we first tested both ordinary least squares regression and robust linear regression. When using ordinary least squares regression, we also tested Cook’s distance filtering to remove outliers at a threshold of Cook’s distance <n/4 (n is the number of datapoints in the regression analysis). A coefficient of variation threshold of 250 was employed as a quality control step in this scaling process. Due to the lower number of false positives in the approach with ordinary least squares regression combined with Cook’s distance filtering, we continued with this approach for all analyses but also report the findings from using the robust linear regression analysis (see Table S3 at https://doi.org/10.6084/m9.figshare.22280632). A detailed workflow for sequin normalization is provided in the vignette for the SIPmg R package (https://github.com/ZielsLab/SIPmg).
In addition to sequin-based normalization, we also explored genome abundance estimation using: (i) unscaled coverage; (ii) relative coverage; and (iii) absolute abundance as per the approach of Greenlon and colleagues (25); and as per the approach of Starr and colleagues (24). Unscaled coverages represented raw average MAG coverage values that were directly used in the estimation of mean weighted BDs and AFE. Relative coverage was estimated as follows: (coverage of an MAG within a fraction)/(sum of coverages of all MAGs within a fraction).

Estimation of atom fraction excess of MAGs

The qSIP model (Equation 1) or ΔBD model (Equation 6) can be used to estimate the AFE of genomes. Briefly, the AFE of organism i can be quantified using the qSIP approach (11):
AFEC,i= MLab,i  MLight,iMHeavymax,i  MLight,i. (1  0.01111233)
(Equation 1-A)
AFEO,i= MLab,i  MLight,iMHeavymax,i  MLight,i. (1  0.002000429)
(Equation 1-B)
where AC,i and AO,i are the estimated AFE with oxygen and carbon as the isotopic substrate, respectively. M Light is the molecular weight of an MAG (g/mol) in the control condition (Equation 2), MLab is the molecular weight of an MAG (g/mol) in the treatment condition (Equation 3), and MHeavymax is the theoretical maximum molecular weight of an MAG (g/mol) due to the maximum labeling by the heavy isotope (Equation 4) in the treatment condition:
MLight= 0.496Gi+307.691
(Equation 2)
MLab= MLight . (WLab WLightWLight+1)
(Equation 3)
MHeavymax= MLight+ 9.974564  0.4987282Gi
(Equation 4)
where Gi is the GC content of the MAG (ranging from 0 to 1). Here, we modified the qSIP model to use the GC content values of MAGs provided from output of CheckM (60), rather than inferring it using an empirical regression (11). WLight and WLab are the mean weighted buoyant densities across control and treatment conditions, respectively.
The weighted average buoyant density (Wij) is then estimated as
Wij = k=1kρjk.  yijkyij
(Equation 5)
where ρjk is the buoyant density of fraction k in replicate j, yijk is the absolute concentration of taxon i in fraction k of replicate j, and yij is the sum total of absolute concentration of taxon i in replicate j. Here, genome abundances were determined using (i) sequin normalization; (ii) relative abundance per coverage and/or reads mapped multiplied by total DNA concentrations; and (iii) relative coverage. Equation 5 above for weighted buoyant density produces a single metric regardless of whether the DNA concentration distribution along the density gradient is monomodal or multimodal.
The estimation of AFE based on the ΔBD model can be represented as
AFEΔBD = WLab WLightImax
(Equation 6)
where Imax is the maximum linear shift in DNA BD (upon 100% labeling), as discussed by Birnie and Rickwood (63). The weighted mean BDs were the same as estimated from Equation 5. This is a variant of ΔBD from the Pepe-Ranney and colleagues’ study (9), in which OTU read counts were interpolated at specific points of the replicate BD gradients to estimate weighted mean BDs. The above models for determining AFE were incorporated into the SIPmg R package for application with SIP metagenomics datasets.

Identifying isotope incorporators using HR-SIP and MW-HR-SIP

To run the high-resolution SIP (HR-SIP) and moving-window HR-SIP (MW-HR-SIP) methods, we used the MAG abundances obtained from the sequin normalization approach. Differential abundances based on absolute abundance for MAGs in the heavy fractions in the treatment conditions were compared to control conditions using HR-SIP and MW-HR-SIP using the HTSSIP R package (31). For HR-SIP, a heavy BD window was set from 1.71 g/mL (as the theoretical peak of E. coli would be at 1.709 g/mL based on a GC content of 0.504) to the maximum buoyant density in every treatment condition. For MW-HR-SIP, the overlapping heavy buoyant density windows chosen were 1.71–1.74 g/mL, 1.72–1.75 g/mL, and 1.73–1.76 g/mL. In all cases, sparsity thresholds between 0% and 30% at 5% intervals were chosen (e.g., sparsity threshold of 25% maintains that MAGs must be present in >25% of fractions in the testing windows). The sparsity threshold with the greatest number of rejected hypotheses was selected for final inference of incorporator identity. The Benjamini–Hochberg method was used to adjust for multiple testing with a threshold of P-value of 0.05 to identify incorporators.

Subsampling of E. coli reads

Reads that mapped to E. coli MAG were extracted from .bam files and subsampled using samtools (v1.7) (htslib 1.7) at 10%, 1%, 0.1%, 0.01%, and 0.001%. New E. coli MAG coverages for each fraction were then calculated (Table S4 at https://doi.org/10.6084/m9.figshare.22280632) and used in SIP analysis to establish limitations that low coverage input may have on the efficiency of bacterial incorporator identification.


The work (proposal: 10.46936/10.25585/60001104) conducted by the U.S. Department of Energy Joint Genome Institute (https://ror.org/04xm1d337), a DOE Office of Science User Facility, is supported by the Office of Science of the U.S. Department of Energy under contract no. DE-AC02-05CH11231. P.S. was supported by a Discovery Grant to R.Z. (RGPIN-2018-04585) from the Natural Sciences and Engineering Research Council of Canada (NSERC). E.E.N., S.J.B., and J.P.R. were supported by a JGI Emerging Technologies Opportunities Program award (DOI: 10.46936/10.25585/60008401) and DOE award SWC1632. Work at Lawrence Livermore National Laboratory was conducted under the auspices of the U.S. DOE under Contract DE-AC52-07NA27344. We thank Tim Mercer at the Garvan Institute of Medical Research for providing sequins.
E.F., R.R.M., R.M.Z., E.E.N., S.J.B., J.P.R., and D.V. conceived the study. D.V. was in charge of overall direction and planning. K.S., A.T., E.E.N., K.B.L., and W.B.K. carried out laboratory work. R.Z., P.S., and D.V. developed the analytical strategies for utilizing internal standards for quantitative SIP metagenomics, and D.V. and P.S. performed the computational analysis. P.S. developed SIPmg R-package. A.C. and R.R. performed MetaHipMer co-assembly. N.V. performed metagenome binning and ANI-AF computation. S.R. contributed to the interpretation of the results. M.K. supported metagenome binning coverage analysis. D.V. and P.S. wrote the manuscript. E.F., R.R.M., and R.M.Z. supervised the study. All authors provided critical feedback and helped shape the research, analysis, and manuscript.


Figure S1 - msystems.01280-22-s0001.eps
A) Average completeness and average purity of MAGs grouped by assembly type (I - intact metagenome assembly with MetaSPAdes, F - separate fractions assembled with metaSPAdes, S - all fractions within each replicate co-assembled with metaSPAdes, M - combined assembly of all fractions using MetaHipMer(v. B) Average coverage across all fractions for each medium and high-quality MAG. Color-coding identifies MAGs found in multiple assembly types (Shared) or uniquely generated in one of the three different assembly types (F - separate fractions assembled with metaSPAdes, S - all fractions within each replicate co-assembled with metaSPAdes, M - combined assembly of all fractions using MetaHipMer). Assemblies of unfractionated DNA (Intact) with MetaSPAdes did not generate unique MAGs.
Figure S2 - msystems.01280-22-s0002.eps
Detecting anomalous samples using pre-centrifugation spike-ins. A) SIP sample displaying the expected spike-in distribution pattern based on relativized absolute coverage along the density gradient. B) An anomalous sample whose spike-in patterns do not match expectations, indicating possible problems in gradient collection and library creation.
Figure S3 - msystems.01280-22-s0003.pdf
Linear regression parameters and Spearman correlations between estimated and expected AFEs obtained using the modified qSIP model (A-H) and the ΔBD method (I-P). (A and I) raw coverage, (B and J) relative coverage, (C and K) multiplying relative abundance with DNA concentration following Greenlon and colleagues (25), (D and L) multiplying relative coverage with DNA concentration following Starr and colleagues (24), (E and M) Sequin approach with ordinary least squares regression without Cook's distance filtering (F and N) Sequin approach with ordinary least squares regression with Cook's distance filtering (G and O) Sequin approach with robust linear regression, and (H and P) Relativizing abundances per fraction (MAG abundance/sum of MAG abundances in each fraction) from sequin approach with robust linear regression. preg and pcor correspond to the P-values for the linear regression and Spearman correlation. The intercepts determined by linear regression were not significantly different from zero (P-value > 0.05) in any method for estimating abundance.
Figure S4 - msystems.01280-22-s0004.pdf
Linear regression parameters and Spearman correlations between estimated and expected AFEs obtained using the qSIP method for subsampled data at mean cumulative coverages of (A) 0.01X, (B) 0.1X, (C) 1X, (D) 10X, (E) 100X, and (F) 1000X. preg and pcor correspond to the P-values for the linear regression and Spearman correlation. The intercepts determined by linear regression were not significantly different from zero (P-value > 0.05) at any level of subsampling.
Figure S5 - msystems.01280-22-s0005.eps
Mean total coverage of MAGs across biological replicates in the unlabeled controls. False positive MAGs are indicated by blue bars (also indicated by arrows). The mean coverage threshold where false positives would be removed (17X) is indicated by a dashed horizontal line. A total of 68 MAGs had mean total coverages greater than this threshold. MAGs lower than this threshold are separated by a dashed vertical line.
Data Set S1 - msystems.01280-22-s0006.xlsx
Internal calibration standards utilized in experimental design. A set of six synthetic DNA fragments (pre) were added to mixtures of DNA from isolates and the complex microbiome to track the ultracentrifugation and fraction collection steps. An additional set of 80 synthetic DNA fragments (post), or sequins, were added to each fraction as an internal standard for subsequent quantitative metagenomic analysis.
Data Set S2 - msystems.01280-22-s0007.xlsx
Metagenome-assembled genomes (MAGs) generated across assembly approaches and associated quality metrics. A total of 2,022 MAGs were generated across all assemblies, of which 248 were high- quality, 447 were medium-quality, and 1,327 were low-quality as defined by the MIMAG reporting standards. Bin identifiers and assembly identifiers are provided, along with CheckM metrics for estimates of completeness and contamination. Cluster representatives are denoted based on single-linkage clustering from average nucleotide identity values of ≥ 96.5 and alignment fractions of ≥ 30%.
Data Set S3 - msystems.01280-22-s0008.xlsx
Metagenome assembly types, metrics, and associated accessions for GOLD and NCBI.
Supplemental Text - msystems.01280-22-s0009.docx
Additional methodological details regarding experimental design, sequencing, assembly, and data analysis.
ASM does not own the copyrights to Supplemental Material that may be linked to, or accessed through, an article. The authors have granted ASM a non-exclusive, world-wide license to publish the Supplemental Material files. Please contact the corresponding author directly for reuse.


[This article was published on 28 June 2023 with incorrect values in the last column of Table 1, inaccurate wording of two sentences and an incorrect citation of supplemental material in the Results section, and inaccuracies in Fig. S3 and S4 in the supplemental material. These items were corrected in the current version, posted on 7 July 2023.]


Falkowski PG, Fenchel T, Delong EF. 2008. The microbial engines that drive Earth’s biogeochemical cycles. Science 320:1034–1039.
Rappé MS, Giovannoni SJ. 2003. The uncultured microbial majority. Annu Rev Microbiol 57:369–394.
Dumont MG, Murrell JC. 2005. Stable isotope probing-linking microbial identity to function. Nat Rev Microbiol 3:499–504.
Radajewski S, Ineson P, Parekh NR, Murrell JC. 2000. Stable-isotope probing as a tool in microbial ecology. Nature 403:646–649.
Coyotzi S, Pratscher J, Murrell JC, Neufeld JD. 2016. Targeted metagenomics of active microbial populations with stable-isotope probing. Curr Opin Biotechnol 41:1–8.
Neufeld JD, Vohra J, Dumont MG, Lueders T, Manefield M, Friedrich MW, Murrell JC. 2007. DNA stable-isotope probing. Nat Protoc 2:860–866.
Neufeld JD, Dumont MG, Vohra J, Murrell JC. 2007. Methodological considerations for the use of stable isotope probing in microbial ecology. Microb Ecol 53:435–442.
Dunford EA, Neufeld JD. 2010. DNA stable-isotope probing (DNA-SIP). J Vis Exp:2027.
Pepe-Ranney C, Campbell AN, Koechli CN, Berthrong S, Buckley DH. 2016. Unearthing the ecology of soil microorganisms using a high resolution DNA-SIP approach to explore cellulose and xylose metabolism in soil. Front Microbiol 7:703.
Barnett SE, Youngblut ND, Buckley DH. 2019. Data analysis for DNA stable isotope probing experiments using multiple window high-resolution SIP. Methods Mol Biol 2046:109–128.
Hungate BA, Mau RL, Schwartz E, Caporaso JG, Dijkstra P, van Gestel N, Koch BJ, Liu CM, McHugh TA, Marks JC, Morrissey EM, Price LB. 2015. Quantitative microbial ecology through stable isotope probing. Appl Environ Microbiol 81:7570–7581.
Buckley DH, Huangyutitham V, Hsu SF, Nelson TA. 2007. Stable isotope probing with 15N achieved by disentangling the effects of genome G+C content and isotope enrichment on DNA density. Appl Environ Microbiol 73:3189–3195.
Long AM, Hou S, Ignacio-Espinoza JC, Fuhrman JA. 2021. Benchmarking microbial growth rate predictions from metagenomes. ISME J 15:183–195.
Purcell AM, Dijkstra P, Finley B, Hayer M, Koch BJ, Mau RL, Morrissey E, Papp K, Schwartz E, Stone BW, Hungate BA. 2020. Quantitative stable Isotope probing with H2 18O to measure taxon-specific microbial growth. Soil Sci Soc Am J 84:1503–1518.
Blazewicz SJ, Hungate BA, Koch BJ, Nuccio EE, Morrissey E, Brodie EL, Schwartz E, Pett-Ridge J, Firestone MK. 2020. Taxon-specific microbial growth and mortality patterns reveal distinct temporal population responses to rewetting in a California grassland soil. ISME J 14:1520–1532.
Kapustina Ž, Medžiūnė J, Alzbutas G, Rokaitis I, Matjošaitis K, Mackevičius G, Žeimytė S, Karpus L, Lubys A. 2021. High-resolution microbiome analysis enabled by linking of 16S rRNA gene sequences with adjacent genomic contexts. Microb Genom 7:000624.
Pinnell LJ, Charles TC, Neufeld JD. 2010. Stable isotope probing and metagenomics, p 97–114. In Stable isotope probing and related technologies.
Ziels RM, Sousa DZ, Stensel HD, Beck DAC. 2018. DNA-SIP based genome-centric metagenomics identifies key long-chain fatty acid-degrading populations in anaerobic digesters with different feeding frequencies. ISME J 12:112–123.
Thomas F, Corre E, Cébron A. 2019. Stable isotope probing and metagenomics highlight the effect of plants on uncultured phenanthrene-degrading bacterial consortium in polluted soil. ISME J 13:1814–1830.
Starr EP, Nuccio EE, Pett-Ridge J, Banfield JF, Firestone MK. 2019. Metatranscriptomic reconstruction reveals RNA viruses with the potential to shape carbon cycling in soil. Proc Natl Acad Sci U S A 116:25900–25908.
Wilhelm RC, DeRito CM, Shapleigh JP, Madsen EL, Buckley DH. 2021. Phenolic acid-degrading Paraburkholderia prime decomposition in forest soil. ISME Commun 1:4.
Sieradzki ET, Morando M, Fuhrman JA. 2021. Metagenomics and quantitative stable isotope probing offer insights into metabolism of polycyclic aromatic hydrocarbon degraders in chronically polluted seawater. mSystems 6:e00245-21.
Yu J, Pavia MJ, Deem LM, Crow SE, Deenik JL, Penton CR. 2020. DNA-stable isotope probing shotgun metagenomics reveals the resilience of active microbial communities to biochar amendment in oxisol soil. Front Microbiol 11:587972.
Starr EP, Shi S, Blazewicz SJ, Koch BJ, Probst AJ, Hungate BA, Pett-Ridge J, Firestone MK, Banfield JF. 2021. Stable-isotope-informed, genome-resolved metagenomics uncovers potential cross-kingdom interactions in rhizosphere soil. mSphere 6:e0008521.
Greenlon A, Sieradzki E, Zablocki O, Koch BJ, Foley MM, Kimbrel JA, Hungate BA, Blazewicz SJ, Nuccio EE, Sun CL, Chew A, Mancilla C-J, Sullivan MB, Firestone M, Pett-Ridge J, Banfield JF. 2022. Quantitative stable-isotope probing (qSIP) with metagenomics links microbial physiology and activity to soil moisture in mediterranean-climate grassland ecosystems. mSystems 7:e0041722.
Nuccio EE, Blazewicz SJ, Lafler M, Campbell AN, Kakouridis A, Kimbrel JA, Wollard J, Vyshenska D, Riley R, Tomatsu A, Hestrin R, Malmstrom RR, Firestone M, Pett-Ridge J. 2022. HT-SIP: a semi-automated stable isotope probing pipeline identifies cross-kingdom interactions in the hyphosphere of arbuscular mycorrhizal fungi. Microbiome 10:199.
Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ. 2017. Microbiome datasets are compositional: and this is not optional. Front Microbiol 8:2224.
McLaren MR, Willis AD, Callahan BJ. 2019. Consistent and correctable bias in metagenomic sequencing experiments. Elife 8:e46923.
Nayfach S, Pollard KS. 2016. Toward accurate and quantitative comparative metagenomics. Cell 166:1103–1116.
Hardwick SA, Chen WY, Wong T, Kanakamedala BS, Deveson IW, Ongley SE, Santini NS, Marcellin E, Smith MA, Nielsen LK, Lovelock CE, Neilan BA, Mercer TR. 2018. Synthetic microbe communities provide internal reference standards for metagenome sequencing and analysis. Nat Commun 9:3096.
Youngblut ND, Barnett SE, Buckley DH. 2018. HTSSIP: an R package for analysis of high throughput sequencing data from nucleic acid stable isotope probing (SIP) experiments. PLoS One 13:e0189616.
Benjamini Y, Yekutieli D. 2005. False discovery rate–adjusted multiple confidence intervals for selected parameters. J Am Stat Assoc 100:71–81.
Hofmeyr S, Egan R, Georganas E, Copeland AC, Riley R, Clum A, Eloe-Fadrosh E, Roux S, Goltsman E, Buluç A, Rokhsar D, Oliker L, Yelick K. 2020. Terabase-scale metagenome coassembly with MetaHipMer. Sci Rep 10:10689.
Kang DD, Li F, Kirton E, Thomas A, Egan R, An H, Wang Z. 2019. MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ 7:e7359.
Bowers RM, Kyrpides NC, Stepanauskas R, Harmon-Smith M, Doud D, Reddy TBK, Schulz F, Jarett J, Rivers AR, Eloe-Fadrosh EA, Tringe SG, Ivanova NN, Copeland A, Clum A, Becraft ED, Malmstrom RR, Birren B, Podar M, Bork P, Weinstock GM, Garrity GM, Dodsworth JA, Yooseph S, Sutton G, Glöckner FO, Gilbert JA, Nelson WC, Hallam SJ, Jungbluth SP, Ettema TJG, Tighe S, Konstantinidis KT, Liu W-T, Baker BJ, Rattei T, Eisen JA, Hedlund B, McMahon KD, Fierer N, Knight R, Finn R, Cochrane G, Karsch-Mizrachi I, Tyson GW, Rinke C, Lapidus A, Meyer F, Yilmaz P, Parks DH, Eren AM, Schriml L, Banfield JF, Hugenholtz P, Woyke T, Genome Standards Consortium. 2017. Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. Nat Biotechnol 35:725–731.
Varghese NJ, Mukherjee S, Ivanova N, Konstantinidis KT, Mavrommatis K, Kyrpides NC, Pati A. 2015. Microbial species delineation using whole genome sequences. Nucleic Acids Res 43:6761–6771.
Koch BJ, McHugh TA, Hayer M, Schwartz E, Blazewicz SJ, Dijkstra P, Gestel N, Marks JC, Mau RL, Morrissey EM, Pett‐Ridge J, Hungate BA. 2018. Estimating taxon‐specific population dynamics in diverse microbial communities. Ecosphere 9:e02090.
Satinsky BM, Gifford SM, Crump BC, Moran MA. 2013. Use of internal standards for quantitative metatranscriptome and metagenome analysis. Methods Enzymol 531:237–250.
Lynch MDJ, Neufeld JD. 2015. Ecology and exploration of the rare biosphere. Nat Rev Microbiol 13:217–229.
Chen Y, Dumont MG, Neufeld JD, Bodrossy L, Stralis-Pavese N, McNamara NP, Ostle N, Briones MJI, Murrell JC. 2008. Revealing the uncultivated majority: combining DNA stable-isotope probing, multiple displacement amplification and metagenomic analyses of uncultivated Methylocystis in acidic peatlands. Environ Microbiol 10:2609–2622.
Sieradzki ET, Koch BJ, Greenlon A, Sachdeva R, Malmstrom RR, Mau RL, Blazewicz SJ, Firestone MK, Hofmockel KS, Schwartz E, Hungate BA, Pett-Ridge J. 2020. Measurement error and resolution in quantitative stable isotope probing: implications for experimental design. mSystems 5:e00151-20.
Sieradzki ET, Greenlon A, Nicolas AM, Firestone MK, Pett-Ridge J, Blazewicz SJ, Banfield JF. 2022. Functional succession of actively growing soil microorganisms during rewetting is shaped by precipitation history. BioRxiv.
Trubl G, Kimbrel JA, Liquet-Gonzalez J, Nuccio EE, Weber PK, Pett-Ridge J, Jansson JK, Waldrop MP, Blazewicz SJ. 2021. Active virus-host interactions at sub-freezing temperatures in Arctic peat soil. Microbiome 9:208.
Starr EP, Shi S, Blazewicz SJ, Probst AJ, Herman DJ, Firestone MK, Banfield JF. 2018. Stable isotope informed genome-resolved metagenomics reveals that Saccharibacteria utilize microbially-processed plant-derived carbon. Microbiome 6:122.
Barnett SE, Buckley DH. 2020. Simulating metagenomic stable isotope probing datasets with MetaSIPSim. BMC Bioinformatics 21:37.
Olm MR, Brown CT, Brooks B, Banfield JF. 2017. DRep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication. ISME J 11:2864–2868.
Sczyrba A, Hofmann P, Belmann P, Koslicki D, Janssen S, Dröge J, Gregor I, Majda S, Fiedler J, Dahms E, Bremges A, Fritz A, Garrido-Oter R, Jørgensen TS, Shapiro N, Blood PD, Gurevich A, Bai Y, Turaev D, DeMaere MZ, Chikhi R, Nagarajan N, Quince C, Meyer F, Balvočiūtė M, Hansen LH, Sørensen SJ, Chia BKH, Denis B, Froula JL, Wang Z, Egan R, Don Kang D, Cook JJ, Deltel C, Beckstette M, Lemaitre C, Peterlongo P, Rizk G, Lavenier D, Wu Y-W, Singer SW, Jain C, Strous M, Klingenberg H, Meinicke P, Barton MD, Lingner T, Lin H-H, Liao Y-C, Silva GGZ, Cuevas DA, Edwards RA, Saha S, Piro VC, Renard BY, Pop M, Klenk H-P, Göker M, Kyrpides NC, Woyke T, Vorholt JA, Schulze-Lefert P, Rubin EM, Darling AE, Rattei T, McHardy AC. 2017. Critical assessment of metagenome Interpretation-a benchmark of metagenomics software. Nat Methods 14:1063–1071.
Meyer F, Fritz A, Deng Z-L, Koslicki D, Lesker TR, Gurevich A, Robertson G, Alser M, Antipov D, Beghini F, Bertrand D, Brito JJ, Brown CT, Buchmann J, Buluç A, Chen B, Chikhi R, Clausen PTLC, Cristian A, Dabrowski PW, Darling AE, Egan R, Eskin E, Georganas E, Goltsman E, Gray MA, Hansen LH, Hofmeyr S, Huang P, Irber L, Jia H, Jørgensen TS, Kieser SD, Klemetsen T, Kola A, Kolmogorov M, Korobeynikov A, Kwan J, LaPierre N, Lemaitre C, Li C, Limasset A, Malcher-Miranda F, Mangul S, Marcelino VR, Marchet C, Marijon P, Meleshko D, Mende DR, Milanese A, Nagarajan N, Nissen J, Nurk S, Oliker L, Paoli L, Peterlongo P, Piro VC, Porter JS, Rasmussen S, Rees ER, Reinert K, Renard B, Robertsen EM, Rosen GL, Ruscheweyh H-J, Sarwal V, Segata N, Seiler E, Shi L, Sun F, Sunagawa S, Sørensen SJ, Thomas A, Tong C, Trajkovski M, Tremblay J, Uritskiy G, Vicedomini R, Wang Z, Wang Z, Wang Z, Warren A, Willassen NP, Yelick K, You R, Zeller G, Zhao Z, Zhu S, Zhu J, Garrido-Oter R, Gastmeier P, Hacquard S, Häußler S, Khaledi A, Maechler F, Mesny F, Radutoiu S, Schulze-Lefert P, Smit N, Strowig T, Bremges A, Sczyrba A, McHardy AC. 2022. Critical assessment of metagenome interpretation: the second round of challenges. Nat Methods 19:429–440.
Quinn TP, Erb I, Gloor G, Notredame C, Richardson MF, Crowley TM. 2019. A field guide for the compositional analysis of any-omics data. Gigascience 8:giz107.
Crossette E, Gumm J, Langenfeld K, Raskin L, Duhaime M, Wigginton K. 2021. Metagenomic quantification of genes with internal standards. mBio 12:e03173-20.
Gifford SM, Sharma S, Rinta-Kanto JM, Moran MA. 2011. Quantitative analysis of a deeply sequenced marine microbial metatranscriptome. ISME J 5:461–472.
Zaramela LS, Tjuanta M, Moyne O, Neal M, Zengler K. 2022. SynDNA—a synthetic DNA spike-in method for absolute quantification of shotgun metagenomic sequencing. mSystems 7:e0044722.
Langenfeld K, Hegarty B, Vidaurri S, Crossette E, Duhaime M, Wigginton K. 2022. A quantitative metagenomic approach to determine population concentrations with examination of quantitative limitations. BioRxiv.
Mondo SJ, Dannebaum RO, Kuo RC, Louie KB, Bewick AJ, LaButti K, Haridas S, Kuo A, Salamov A, Ahrendt SR, Lau R, Bowen BP, Lipzen A, Sullivan W, Andreopoulos BB, Clum A, Lindquist E, Daum C, Northen TR, Kunde-Ramamoorthy G, Schmitz RJ, Gryganskyi A, Culley D, Magnuson J, James TY, O’Malley MA, Stajich JE, Spatafora JW, Visel A, Grigoriev IV. 2017. Widespread adenine N6-methylation of active genes in fungi. Nat Genet 49:964–968.
Farthing DE, Buxbaum NP, Lucas PJ, Maglakelidze N, Oliver B, Wang J, Hu K, Castro E, Bare CV, Gress RE. 2017. Comparing DNA enrichment of proliferating cells following administration of different stable isotopes of heavy water. Sci Rep 7:4043.
Nurk S, Meleshko D, Korobeynikov A, Pevzner PA. 2017. metaSPAdes: a new versatile metagenomic assembler. Genome Res 27:824–834.
Mikheenko A, Saveliev V, Gurevich A. 2016. MetaQUAST: evaluation of metagenome assemblies. Bioinformatics 32:1088–1090.
Clum A, Huntemann M, Bushnell B, Foster B, Foster B, Roux S, Hajek PP, Varghese N, Mukherjee S, Reddy TBK, Daum C, Yoshinaga Y, O’Malley R, Seshadri R, Kyrpides NC, Eloe-Fadrosh EA, Chen I-M, Copeland A, Ivanova NN. 2021. DOE JGI metagenome workflow. mSystems 6:e00804-20.
Chaumeil P-A, Mussig AJ, Hugenholtz P, Parks DH. 2022. GTDB-Tk v2: memory friendly classification with the genome taxonomy database. Bioinformatics 38:5315–5316.
Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. 2015. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res 25:1043–1055.
Mikheenko A, Prjibelski A, Saveliev V, Antipov D, Gurevich A. 2018. Versatile genome assembly evaluation with QUAST-LG. Bioinformatics 34:i142–i150.
Chen I-MA, Chu K, Palaniappan K, Ratner A, Huang J, Huntemann M, Hajek P, Ritter SJ, Webb C, Wu D, Varghese NJ, Reddy TBK, Mukherjee S, Ovchinnikova G, Nolan M, Seshadri R, Roux S, Visel A, Woyke T, Eloe-Fadrosh EA, Kyrpides NC, Ivanova NN. 2023. The IMG/M data management and analysis system v.7: content updates and new features. Nucleic Acids Res 51:D723–D732.
Summers WC. 1979. Centrifugal separations in molecular and cell biology. Yale J Biol Med 52:491–492.

Information & Contributors


Published In

cover image mSystems
Volume 8Number 431 August 2023
eLocator: e01280-22
Editor: Rachel Poretsky, University of Illinois at Chicago, Chicago, Illinois, USA
PubMed: 37377419


Received: 19 December 2022
Accepted: 19 April 2023
Published online: 28 June 2023


  1. stable isotope probing
  2. metagenomics
  3. DNA-SIP
  4. co-assembly
  5. internal standards
  6. spike-ins

Data Availability

Raw metagenome sequencing reads have been deposited under BioProject Accession PRJNA878529. The MetaHipMer combined assembly and annotated data is available through IMG/M under taxon identifier 3300045762. Single-fraction and combined per-sample assemblies, along with all MAGs and input files for qSIP analysis are available via https://portal.nersc.gov/dna/microbial/prokpubs/DVyshenska2022/. A full list of available data and associated NCBI accession numbers are available in Data Set S3.
The code for the SIPmg R package is available for download, along with a vignette describing all functions, at: https://github.com/ZielsLab/SIPmg. The SIPmg package includes functions to calculate global scaling factors for genomes based on regression of sequin coverage versus concentration using either ordinary least squares linear regression or robust linear regression. The package can thereafter estimate AFE using either qSIP model or ΔBD method. The package also outputs both FDR adjusted and Bonferroni adjusted bootstrapped AFE confidence intervals for MAGs. The package can also perform HR-SIP and MW-HR-SIP which were built using the HTS-SIP R package.



Dariia Vyshenska
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contributions: Conceptualization, Data curation, Formal analysis, Methodology, Project administration, Validation, Visualization, Writing – original draft, and Writing – review and editing.
Pranav Sampara
Department of Civil Engineering, The University of British Columbia, Vancouver, British Columbia, Canada
Author Contributions: Formal analysis, Software, Validation, Visualization, Writing – original draft, and Writing – review and editing.
Kanwar Singh
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contribution: Methodology.
Andy Tomatsu
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contribution: Methodology.
W. Berkeley Kauffman
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contribution: Methodology.
Erin E. Nuccio
Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, California, USA
Author Contributions: Conceptualization and Methodology.
Steven J. Blazewicz
Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, California, USA
Author Contribution: Conceptualization.
Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, California, USA
Life & Environmental Sciences Department, University of California Merced, Merced, California, USA
Author Contributions: Conceptualization, Supervision, and Writing – review and editing.
Katherine B. Louie
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contribution: Methodology.
Neha Varghese
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contribution: Formal analysis.
Matthew Kellom
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contribution: Formal analysis.
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contribution: Formal analysis.
Robert Riley
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contribution: Formal analysis.
Simon Roux
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contribution: Validation.
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contributions: Conceptualization, Supervision, and Writing – review and editing.
Department of Civil Engineering, The University of British Columbia, Vancouver, British Columbia, Canada
Author Contributions: Conceptualization, Funding acquisition, Software, Supervision, and Writing – review and editing.
DOE Joint Genome Institute, Lawrence Berkeley National Laboratory, Berkeley, California, USA
Author Contributions: Conceptualization, Supervision, and Writing – review and editing.


Rachel Poretsky
University of Illinois at Chicago, Chicago, Illinois, USA


Dariia Vyshenska and Pranav Sampara contributed equally to this work. Author order was determined by the corresponding authors after negotiation.
Ryan M. Ziels and Rex R. Malmstrom contributed equally to this work.
The authors declare no conflict of interest.

Metrics & Citations



  • For recently published articles, the TOTAL download count will appear as zero until a new month starts.
  • There is a 3- to 4-day delay in article usage, so article usage will not appear immediately after publication.
  • Citation counts come from the Crossref Cited by service.


If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. For an editable text file, please select Medlars format which will download as a .txt file. Simply select your manager software from the list below and click Download.

View Options

Figures and Media






Share the article link

Share with email

Email a colleague

Share on social media

American Society for Microbiology ("ASM") is committed to maintaining your confidence and trust with respect to the information we collect from you on websites owned and operated by ASM ("ASM Web Sites") and other sources. This Privacy Policy sets forth the information we collect about you, how we use this information and the choices you have about how we use such information.
FIND OUT MORE about the privacy policy