Metagenomic Information Recovery from Human Stool Samples Is Influenced by Sequencing Depth and Profiling Method

Sequencing of the 16S rRNA gene (16S) has long been a go-to method for microbiome characterization due to its accessibility and lower cost compared to shotgun metagenomic sequencing (SMS). However, 16S sequencing rarely provides species-level resolution and cannot provide direct assessment of other taxa (e.g., viruses and fungi) or functional gene content. Shallow shotgun metagenomic sequencing (SSMS) has emerged as an approach to bridge the gap between 16S sequencing and deep metagenomic sequencing. SSMS is cost-competitive with 16S sequencing, while also providing species-level resolution and functional gene content insights. In the present study, we evaluated the effects of sequencing depth on marker gene-mapping- and alignment-based annotation of bacteria in healthy human stool samples. The number of identified taxa decreased with lower sequencing depths, particularly with the marker gene-mapping-based approach. Other annotations, including viruses and pathways, also showed a depth-dependent effect on feature recovery. These results refine the understanding of the suitability and shortcomings of SSMS, as well as annotation tools for metagenomic analyses in human stool samples. Results may also translate to other sample types and may open the opportunity to explore the effect of sequencing depth and annotation method.


Introduction
Sequencing-based studies are commonly employed to identify patterns in the microbiome and cause-effect relationships of microorganisms with their host and environment [1,2]. Although both 16S rRNA gene (16S) sequencing and shotgun metagenomic sequencing (SMS) can be utilized for this purpose, traditional SMS is often more cost-prohibitive than 16S sequencing and typically requires substantially greater computational resources to process [3]. For these reasons, 16S sequencing is the effect of SMS depth on pathway identification and viral communities, using healthy human stool samples as proxies of the gut microbiome.

Samples Processing for 16S rRNA Sequencing
Ten human stool samples were purchased from Lee Biosolutions (Maryland Heights, MO, USA) and stored at −20 • C until further processing. ATCC ® mock community MSA1001, comprised of DNA from 10 organisms in staggered concentrations (ranging in relative abundances from 44.78 to 0.04%) was processed in parallel as a control. A no template PCR sample was included as a negative control. For 16S rRNA gene sequencing of the ten stool samples, the V4 hypervariable region was amplified by PCR and sequenced on the MiSeq platform (Illumina; San Diego, CA, USA) using the 2 × 250 bp paired-end protocol following manufacturer's instructions (Illumina; San Diego, CA, USA). Amplification of the V4 region of the 16S gene was selected in the present study, because it is well-standardized, commonly utilized, and has been shown to capture greater diversity than other variable regions in stool samples [21]. The primers used for amplification contain adapters for sequencing in the Illumina MiSeq platform, which also contains single-end barcodes, allowing pooling and direct sequencing of PCR products [22]. The read pairs were demultiplexed followed by merging using USEARCH v7.0.1090 [23], allowing zero mismatches and a minimum overlap of 50 bases. Merged reads were trimmed using Trimmomatic (version 0.39) to perform adapter trimming, quality trimming, and filtering (SLIDINGWINDOW:4:20 MINLEN:75) at the first base with Q5 [24]. In addition, a quality filter was applied to the resulting merged reads and reads containing above 0.05 expected errors were discarded. The 16S rRNA gene sequences were clustered into operational taxonomic units (OTUs) at a similarity cutoff value of 97% using the UPARSE algorithm [25]. OTUs were mapped to an optimized version of the SILVA (v132) database containing only the 16S v4 region from 6,087,080 rRNA gene sequences to determine taxonomies. Abundances were recovered by mapping the demultiplexed reads to the UPARSE OTUs. All downstream analyses including α-diversity, β-diversity, and taxonomic analyses were performed using in-house data visualization tools. For comparison with SMS data, OTUs were collapsed to genus level.

Shotgun Metagenomic Library Preparation, Sequencing, and Data Quality
Genomic DNA from the ten stool samples described above, as well as ATCC ® mock community MSA1001, were prepared for sequencing using the Nextera DNA Flex library preparation kit (Illumina; San Diego, CA, USA) with Nextera Index Kit (Illumina; San Diego, CA, USA). Library size estimation and quantification were determined with the fragment analyzer (Advanced Analytical Technologies, Inc., Ankeny, IA, USA) electrophoresis system. Sequencing was performed on a NextSeq using the 2 × 150 bp paired-end protocol following manufacturer's instructions (Illumina; San Diego, CA, USA). After sequencing, Trimmomatic (version 0.39) was used as described above [24]. After quality trimming and filtering, host filtering was performed by aligning reads to a human reference genome (GRCh38; GCA_000001405.15) that was pulled from Bowtie2 (version 2.3.4.1) with default parameters [26]. Unmapped reads were retained for downstream analysis, including taxonomic, pathway, and viral profiling.

Taxonomic and Pathway Annotations
After subsampling the reads, bacterial taxonomic annotations were performed using MetaPhlAn2 (v2.7.7) [13], which uses clade-specific markers that provide bacterial, archaeal, viral, and eukaryotic quantification at the species level using default parameters for reported relative abundances. This version of the MetaPhlAn database was developed using 16,904 reference genomes. Since the bacterial component was the most abundant, only bacteria were considered, and data were re-normalized prior to downstream analysis. DNA sequences were also aligned to the Venti database (originally developed at CoreBiome), a curated proprietary database containing 19,840 complete and/or high-quality draft bacterial genomes from RefSeq. Alignments were made at 97% identity against all reference genomes. Every input sequence was compared to every reference sequence in the Venti database using fully gapped alignment with the BURST mapping algorithm [14]. Ties were broken by minimizing the overall number of unique genome hits. For taxonomy assignment, each input sequence was assigned the lowest common ancestor that was consistent across at least 80% of all reference sequences tied for best hit. The number of counts for each genome was normalized to genome length, and genomes accounting for less than one millionth of all species-level hits and those with less than 0.01% of their unique genome regions covered (and <1% of the whole genome) were discarded. Samples with fewer than 10,000 sequences were also discarded. Count data were converted to relative abundance for each sample. Both MetaPhlAn2 and BURST data were further filtered to remove taxa present in relative abundances <0.01% using the script filter_otus_from_otu_table.py in the Quantitative Insights into Microbial Ecology (QIIME) platform [27]. The normalized and filtered  [28]. VirMap merges nucleotide and protein information to assign taxonomy independently from genome coverage or read overlap. Relative abundances of selected viruses were acquired by dividing the number of counts of the virus of interest per sequencing depth divided by the genome size in bp.
Pathway profiling was performed using HUMAnN2 v0.11.1 with default parameters [29]. Pathway searches were performed against UniProt and MetaCyc50 [30,31]. Unclassified and unmapped categories were removed, and data were re-normalized for further analyses.

Statistical Analyses
Statistical analyses were performed in RStudio version 1.1.456. For α-diversity comparisons, the Friedman Rank test was used to determine statistical differences among sequencing depths, as repeated measure was considered, and the Nemenyi post hoc test was further applied for sequencing depth comparisons. Permutational multivariate analysis of variance (PERMANOVA) was performed using proprietary tools to visualize differences based on weighted Bray-Curtis distances. Distance values were then visualized using Principal Coordinate Analysis (PCoA) plots. Nonparametric permutation tests were employed to test for statistical significance amongst sequencing depths and taxonomy and pathways relative abundances. Briefly, comparisons of genera or species across the various sequencing depths were performed using Kruskal-Wallis in QIIME using the group_significance_analyses.py script, with default parameters. The group_significance.py script compares each feature (e.g., OTU, species, gene family) to see if it is differentially abundant based on the groupings of interest, as defined in the metadata file. The script will group together samples that have the same value in the mapping file under the header passed with the -c option. The script output includes unadjusted p-values, and false discovery rate (FDR) values. Correlation plots were generated using the function ggscatter () in RStudio with the addition of correlation coefficient and p-value to indicate significance.

Data Availability
The sequences generated and analyzed during the current study are available in the NCBI repository under accession number PRJNA676159.  [13], as well as BURST to assess the potential effects of marker gene-mappingand alignment-based approaches on bacterial taxonomic classification, respectively [14]. In addition, we performed 16S sequencing for comparison at the genus level. Prior to further analyses, data were filtered to remove bacterial genera and species with relative abundances < 0.01% to decrease potential false positives. This was performed separately for both stool and mock community samples. α-diversity at the genus level was determined at the various mentioned sequencing depths ( Figure 1).  [13], as well as BURST to assess the potential effects of marker gene-mappingand alignment-based approaches on bacterial taxonomic classification, respectively [14]. In addition, we performed 16S sequencing for comparison at the genus level. Prior to further analyses, data were filtered to remove bacterial genera and species with relative abundances < 0.01% to decrease potential false positives. This was performed separately for both stool and mock community samples. αdiversity at the genus level was determined at the various mentioned sequencing depths ( Figure 1). Results show a slightly lower number of observed genera using SMS and a marker genemapping-based approach at 5 Gb compared to 16S sequencing ( Figure 1A) (Supplementary Dataset 1), although this difference was not statistically significant (p-value = 0.96) (Supplementary Dataset 2). Observed genera significantly decreased at 0.5 Gb compared to 5 Gb (p-value = 0.020) when data were annotated using the marker gene-mapping-based approach (Supplementary Dataset 2). Interestingly, higher numbers of observed genera were noted at all sequencing depths when using the alignment approach compared to 16S sequencing at similar depths ( Figure 1A) (Supplementary Dataset 1). Differences between the alignment vs. the 16S methods were determined using a Friedman rank sum test and were statistically significant at 5 Gb (p-value = 1.40 × 10 −7 ), 3 Gb (p-value = 4.20 × 10 −7 ), 1 Gb (p-value = 8.60 × 10 −4 ), and 0.75 Gb (p-value = 2.65 × 10 −3 ). Observed genera significantly decreased at 0.25 Gb (p-value = 3.79 × 10 −3 ) compared to 5 Gb when data were annotated using the alignment-based approach. Shannon diversity at the genus level did not significantly decrease at any of the sequencing depths tested when using the marker gene-mapping-based-and alignment-based approaches ( Figure 1A and Supplementary Dataset 1). Results show a slightly lower number of observed genera using SMS and a marker gene-mapping-based approach at 5 Gb compared to 16S sequencing ( Figure 1A) (Supplementary Dataset 1), although this difference was not statistically significant (p-value = 0.96) (Supplementary Dataset 2). Observed genera significantly decreased at 0.5 Gb compared to 5 Gb (p-value = 0.020) when data were annotated using the marker gene-mapping-based approach (Supplementary Dataset 2). Interestingly, higher numbers of observed genera were noted at all sequencing depths when using the alignment approach compared to 16S sequencing at similar depths ( Figure 1A) (Supplementary Dataset 1). Differences between the alignment vs. the 16S methods were determined using a Friedman rank sum test and were statistically significant at 5 Gb (p-value = 1.40 × 10 −7 ), 3 Gb (p-value = 4.20 × 10 −7 ), 1 Gb (p-value = 8.60 × 10 −4 ), and 0.75 Gb (p-value = 2.65 × 10 −3 ). Observed genera significantly decreased at 0.25 Gb (p-value = 3.79 × 10 −3 ) compared to 5 Gb when data were annotated using the alignment-based approach. Shannon diversity at the genus level did not significantly decrease at any of the sequencing depths tested when using the marker gene-mapping-based-and alignment-based approaches ( Figure 1A and Supplementary Dataset 1).

Bacterial
α-diversity analyses also included observed species and Shannon diversity at the species level using Friedman rank sum test ( Figure 1B). While the number of observed species was not significantly different between the marker gene-mapping-based approach (5 Gb) and 16S sequencing (p-value = 0.99), the number of observed species identified at 5 Gb (p-value = 1.10 × 10 −8 ), 3 Gb (p-value = 2.00 × 10 −6 ), 1 Gb (p-value = 7.00 × 10 −4 ), and 0.75 Gb (p-value = 0.002) when using the alignment-based approach was significantly higher than those observed when using 16S sequencing. Notably, most of the 16S sequences were unable to be classified at the species level (Supplementary Dataset 3). The number of observed species significantly decreased at 0.5 Gb when using the marker gene-mapping-(p-value = 0.009) and alignment-based approaches (p-value = 0.012) compared to 5 Gb (Supplementary Dataset 2).
β-diversity (Bray-Curtis weighted dissimilarity distances) of the stool samples at the genus level showed separation of samples based on method (marker gene-mapping vs. alignment vs. 16S sequencing), with the 16S sequencing data clustering closer to the alignment-based data (PERMANOVA; p-value = 0.001; R-Squared = 0.338; F-Statistic = 37.5) (Figure 2A). PCoA plots of the marker gene-mapping-based Bray-Curtis distances at the genus level showed clustering by subject rather than sequencing depth (PERMANOVA; p-value = 1; R-Squared = 0.00814; F-Statistic = 0.0862) ( Figure 2B). Interestingly, subject-specific effects were noted, as marker gene-mapping-based data from some subjects at 0.1 Gb did not cluster closely to their counterpart profiles generated at greater sequencing depths ( Figure 2B). Similarly, PCoA plots of the alignment-based Bray-Curtis distances at the genus level showed clustering by subject (PERMANOVA; p-value = 0.155; R-Squared = 3.55 × 10 −6 ; F-Statistic α-diversity analyses also included observed species and Shannon diversity at the species level using Friedman rank sum test ( Figure 1B). While the number of observed species was not significantly different between the marker gene-mapping-based approach (5 Gb) and 16S sequencing (p-value = 0.99), the number of observed species identified at 5 Gb (p-value = 1.10 × 10 −8 ), 3 Gb (pvalue = 2.00 × 10 −6 ), 1 Gb (p-value = 7.00 × 10 −4 ), and 0.75 Gb (p-value = 0.002) when using the alignment-based approach was significantly higher than those observed when using 16S sequencing. Notably, most of the 16S sequences were unable to be classified at the species level (Supplementary Dataset 3). The number of observed species significantly decreased at 0.5 Gb when using the marker gene-mapping-(p-value = 0.009) and alignment-based approaches (p-value = 0.012) compared to 5 Gb (Supplementary Dataset 2).
β-diversity (Bray-Curtis weighted dissimilarity distances) of the stool samples at the genus level showed separation of samples based on method (marker gene-mapping vs. alignment vs. 16S sequencing), with the 16S sequencing data clustering closer to the alignment-based data (PERMANOVA; p-value = 0.001; R-Squared = 0.338; F-Statistic = 37.5) (Figure 2A). PCoA plots of the marker gene-mapping-based Bray-Curtis distances at the genus level showed clustering by subject rather than sequencing depth (PERMANOVA; p-value = 1; R-Squared = 0.00814; F-Statistic = 0.0862) ( Figure 2B). Interestingly, subject-specific effects were noted, as marker gene-mapping-based data from some subjects at 0.1 Gb did not cluster closely to their counterpart profiles generated at greater sequencing depths ( Figure 2B). Similarly, PCoA plots of the alignment-based Bray-Curtis distances at the genus level showed clustering by subject (PERMANOVA; p-value = 0.155; R-Squared = 3.55 × 10 −6 ; F-Statistic = 3.73 × 10 −5 ) ( Figure 2C). Bray-Curtis weighted dissimilarity distances of the stool samples at the species level showed separation of the data based on the method used (marker gene mapping vs. alignment) (PERMANOVA; p-value = 0.001; R-Squared = 0.62; F-Statistic = 226) ( Figure  2D). Plotting of the marker gene-mapping-(PERMANOVA; p-value = 1; R-Squared = 0.0148; F-Statistic = 0.157) ( Figure 2E) and alignment-based data (PERMANOVA; p-value = 1; R-Squared = 1.73 × 10 −5 ; F-Statistic = 0.000182) ( Figure 2F) at the species level showed clustering by subject.  (B) SMS data annotated using the marker gene-mapping method at the genus level. (C) SMS data annotated using the alignment method at the genus level. (D) SMS data annotated using the marker gene mapping and alignment methods at the species level. (E) SMS data annotated using the marker gene-mapping method at the species level. (F) SMS data annotated using the alignment method at the species level. Clustering of data by subjects is highlighted using circles or ellipses.

Identification of Certain Bacterial Genera and Species Is Affected by SMS Sequencing Depth When Using a Marker Gene-Mapping-Based Approach, but Not an Alignment-Based Approach
The effects of marker gene-mapping-and alignment-based tools at various SMS sequencing depths, as well as 16S sequencing, to the relative abundances of bacterial genera were determined (Supplementary Dataset 3 and Supplementary Dataset 4). The average relative abundances of common gut commensals typically varying in abundance in human stool, including Bifidobacterium ( Figure 3A), Alistipes ( Figure 3B), Roseburia ( Figure 3C), and Lactobacillus ( Figure 3D) were plotted. Relative abundances of these commensals were relatively consistent across sequencing depths, with the exception of Lactobacillus, when using the marker gene-mapping-versus alignment-based approaches (Supplementary Dataset 3 and Supplementary Dataset 4). Lactobacillus, which is present at relatively low abundance, was rarely detectable when using the marker gene-mapping-based approach at a sequencing depth < 0.25 Gb ( Figure 3D)

Identification of Certain Bacterial Genera and Species Is Affected by SMS Sequencing Depth When Using a Marker Gene-Mapping-Based Approach, but not an Alignment-Based Approach
The effects of marker gene-mapping-and alignment-based tools at various SMS sequencing depths, as well as 16S sequencing, to the relative abundances of bacterial genera were determined (Supplementary Dataset 3 and Supplementary Dataset 4). The average relative abundances of common gut commensals typically varying in abundance in human stool, including Bifidobacterium ( Figure 3A), Alistipes ( Figure 3B), Roseburia ( Figure 3C), and Lactobacillus ( Figure 3D) were plotted. Relative abundances of these commensals were relatively consistent across sequencing depths, with the exception of Lactobacillus, when using the marker gene-mapping-versus alignment-based approaches (Supplementary Dataset 3 and Supplementary Dataset 4). Lactobacillus, which is present at relatively low abundance, was rarely detectable when using the marker gene-mapping-based approach at a sequencing depth < 0.25 Gb ( Figure 3D  Correlations of all bacterial genera across the various SMS depths using the marker genemapping-and alignment-based approaches, as well as 16S, were performed to demonstrate the extent Correlations of all bacterial genera across the various SMS depths using the marker gene-mappingand alignment-based approaches, as well as 16S, were performed to demonstrate the extent of signal variation relative to depth and annotation approach. Data obtained at 5 Gb from both the marker gene-mapping-and alignment-based approaches showed strong, statistically significant correlations at 3 Gb ( Figure 4A), 1 Gb ( Figure 4B), 0.75 Gb (Figure 4C), 0.5 Gb (Figure 4D), and 0.25 Gb ( Figure 4E) (R =1; p-value < 2.2 × 10 −16 ). Strong, statically significant correlations were also noted with the alignment-based approach at 0.1 Gb ( Figure 4F) (R = 1; p-value < 2.2 × 10 −16 ), but these were slightly decreased when using the marker gene-mapping-based approach (R = 0.98; p-value < 2.2 × 10 −16 ). Moderate (R = 0.67; p-value < 2.2 × 10 −16 ) and weak (R = 0.3; p-value < 2.2 × 10 −16 ) correlations were noted between the alignment-and marker gene-mapping-based information, with 16S sequencing, respectively ( Figure 4G).

Assessing the Accuracy of Marker Gene-Mapping-and Alignment-Based Approaches Using Reference Materials
To assess the accuracy of both the marker gene-mapping-and alignment-based methods, mock community data (ATCC ® MSA1001) were annotated at the various sequencing depths tested. At the genus level, both annotation methods identified all 10 expected taxa. While no false positives were identified with the marker gene-mapping-based approach, the alignment-based approach incorrectly identified Shigella, which was not present in the mock community. The marker gene-mapping-based approach identified all 10 bacterial genera at 5 Gb in relative abundances similar to those expected (Supplementary Dataset 7 and Supplementary Figure S1A). However, Bifidobacterium, Deinococcus, and Enterococcus, expected to be present in relative abundances of 0.04%, were not identified at a sequencing depth < 0.5 Gb. The alignment-based approach underestimated the expected abundance of Escherichia, but the remaining bacterial genera were identified at all sequencing depths in relative abundances compared to the expected relative abundances (Supplementary Dataset 7 and Supplementary Figure S1B).
At the species level, all 10 of the expected taxa were identified when using the marker gene-mapping-based method at 5 Gb. The Deinococcus radiodurans signal was lost at a sequencing depth of <0.5 Gb, while Bifidobacterium adolescentis and Enterococcus faecalis signals were lost at sequencing depths < 0.25 Gb when using the marker gene-mapping-based approach (Supplementary Figure S1C and Supplementary Dataset 7). The alignment-based approach captured all 10 of the expected taxa across the various sequencing depths but underestimated the relative abundance of B. adolescentis and Escherichia coli (Supplementary Figure S1D). Interestingly, approximately 55.4% of the sequences in the mock community were not classified (Unclassified), while 2.5% and 0.1% of the sequences were attributed to Bacillus cereus and Staphylococcus sp. HMSC055B03, respectively (Supplementary Dataset 7).

Sequencing Depth Affects Virus Identification from Metagenomic Samples and Is Subject-Specific
SMS data at 5 Gb, 3 Gb, 1 Gb, 0.75 Gb, 0.5 Gb, 0.25 Gb, and 0.1 Gb were processed for virus taxonomic characterization using VirMap [28]. The total number of identified viruses decreased significantly at a sequencing depth below 0.5 Gb (p-value = 9.0 × 10 −4 ) ( Figure 5A). The stochasticity of sampling sequences with low counts was evident with the virome dataset analyzed. This limited the ability to fully understand the effect of SMS depth relative to viral detection; however, data could still be used as a proof-of-concept. For this, virus abundances from individual, healthy subjects were normalized by dividing the total counts by the virus genome size (bp). Percentages of normalized counts were then plotted. Values for Arthrobacter phage Mendel (taxid 2484218) ( Figure 5B Figure 5F) were plotted. Although not all the viruses shown in Figure 5B-F were identified in all the subjects, these results suggest that 5 Gb recovered more viruses than lower sequencing depths ( Figure 5A), but the implications of these findings should be considered in a subject-to-subject and virus-per-virus scenario (Supplementary Dataset 8).

Pathway Relative Abundances Are Dependent on Sequencing Depth and Are Category-Specific
Pathway information was obtained using HUMAnN2, with UniProt and MetaCyc databases [29]. HUMAnN2 provides a relatively fast and accurate species-resolved functional profiles by aligning the reads to species pangenomes, performing translated search on unclassified reads, and quantifying gene families and pathways. The total number of pathways identified using UniProt significantly decreased at a sequencing depth below 0.75 Gb (p-value = 0.019) ( Figure 6A). The pathways identified using UniProt are shown in Supplementary Dataset 9. The relative abundances of a number of pathways identified using UniProt were not significantly different across the various sequencing depths (Supplementary Dataset 10). Similarly, the total number of pathways identified using MetaCyc were significantly reduced at a sequencing depth below 0.75 Gb (p-value = 0.043) ( Figure 6B). All MetaCyc pathways identified are shown in Supplementary Dataset 11. The number of recovered MetaCyc pathways were significantly affected by sequencing depth (Supplementary Dataset 12). For instance, when looking at the MetaCyc pathway output, L-glutamine biosynthesis pathway was not recovered at a sequencing depth < 0.5 Gb (FDR-adjusted p-value = 9.38 × 10 -4 ). Other categories, such as the superpathway (which usually have an additional preceding class within the pathway ontology in order to define their biological role) of menaquinol-6 biosynthesis (FDRadjusted p-value = 0.020) and the superpathway of demethylmenaquinol-6 biosynthesis (FDR-

Pathway Relative Abundances Are Dependent on Sequencing Depth and Are Category-Specific
Pathway information was obtained using HUMAnN2, with UniProt and MetaCyc databases [29]. HUMAnN2 provides a relatively fast and accurate species-resolved functional profiles by aligning the reads to species pangenomes, performing translated search on unclassified reads, and quantifying gene families and pathways. The total number of pathways identified using UniProt significantly decreased at a sequencing depth below 0.75 Gb (p-value = 0.019) ( Figure 6A). The pathways identified using UniProt are shown in Supplementary Dataset 9. The relative abundances of a number of pathways identified using UniProt were not significantly different across the various sequencing depths (Supplementary Dataset 10). Similarly, the total number of pathways identified using MetaCyc were significantly reduced at a sequencing depth below 0.75 Gb (p-value = 0.043) ( Figure 6B). All MetaCyc pathways identified are shown in Supplementary Dataset 11. The number of recovered MetaCyc pathways were significantly affected by sequencing depth (Supplementary Dataset 12). For instance, when looking at the MetaCyc pathway output, L-glutamine biosynthesis pathway was not recovered at a sequencing depth < 0.5 Gb (FDR-adjusted p-value = 9.38 × 10 −4 ). Other categories, such as the superpathway (which usually have an additional preceding class within the pathway ontology in order to define their biological role) of menaquinol-6 biosynthesis (FDR-adjusted p-value = 0.020) and the superpathway of demethylmenaquinol-6 biosynthesis (FDR-adjusted p-value = 0.020) were not identified at a sequencing depth < 1 Gb (Supplementary Dataset 12).

Discussion
The present study aimed to augment current understanding of the effect of both sequencing depth and annotation methods using human stool samples and a mock community with a goal of helping to inform decision making in terms of sequencing method, depth, and annotation tool, particularly for pilot studies and/or studies with high volume samples. Marker gene-mapping-and alignment-based methods for SMS were applied to determine bacterial taxonomic classification and compared with 16S rRNA gene (V4) genus level profiles. Bacterial taxonomic comparison between SMS and 16S at the genus level and at various sequencing depths has not been extensively evaluated. In addition, bacterial species level classification at the various SMS depths was determined using both marker gene-mapping-and alignment-based tools, which have not extensively been evaluated in parallel. The present study also evaluated the effect of SMS depth on the recovery of functional pathways and viral taxa, two categories of microbiome data that are frequently overlooked. Results of the present study build upon previous studies that have sought to assess the effect of sequencing depth with regard to various body sites using k-mer and alignment methods [32]; stool samples and reference materials using marker, k-mer, and alignment methods [11]; stool and food samples, as well as reference materials using k-mer-based methods [33].
Achieving accurate taxonomic and functional resolution is increasingly important when aiming to identify patterns in microbiome studies. However, this is often not achieved when sequencing one or two 16S variable regions, as with amplicon sequencing. In the present study, α-diversity results at the genus level showed lower numbers of genera when using SMS data annotated using the marker gene-mapping-based approach compared to 16S sequencing. This is, in part, a reflection of the availability of a large number of unclassified OTUs at the genus level in 16S databases, as well as the ability of amplicon-based sequencing to capture low-abundance taxa. In addition, that a number of observed genera were not identified using the marker gene-mapping-or the alignment-based approach indicates that reference databases will continue to benefit from the sequencing of additional novel bacterial. The ability to detect low-abundance bacterial genera and species being impaired at shallow sequencing depths when using the marker gene-mapping-based approach suggests that the combination of shallow sequencing and marker gene-mapping-base annotation may not be suitable for studies seeking to characterize low-abundance bacterial taxa, including those that are transient members of the gut microbiome. Notably, the number of bacterial genera and species detected were

Discussion
The present study aimed to augment current understanding of the effect of both sequencing depth and annotation methods using human stool samples and a mock community with a goal of helping to inform decision making in terms of sequencing method, depth, and annotation tool, particularly for pilot studies and/or studies with high volume samples. Marker gene-mapping-and alignment-based methods for SMS were applied to determine bacterial taxonomic classification and compared with 16S rRNA gene (V4) genus level profiles. Bacterial taxonomic comparison between SMS and 16S at the genus level and at various sequencing depths has not been extensively evaluated. In addition, bacterial species level classification at the various SMS depths was determined using both marker gene-mapping-and alignment-based tools, which have not extensively been evaluated in parallel. The present study also evaluated the effect of SMS depth on the recovery of functional pathways and viral taxa, two categories of microbiome data that are frequently overlooked. Results of the present study build upon previous studies that have sought to assess the effect of sequencing depth with regard to various body sites using k-mer and alignment methods [32]; stool samples and reference materials using marker, k-mer, and alignment methods [11]; stool and food samples, as well as reference materials using k-mer-based methods [33].
Achieving accurate taxonomic and functional resolution is increasingly important when aiming to identify patterns in microbiome studies. However, this is often not achieved when sequencing one or two 16S variable regions, as with amplicon sequencing. In the present study, α-diversity results at the genus level showed lower numbers of genera when using SMS data annotated using the marker gene-mapping-based approach compared to 16S sequencing. This is, in part, a reflection of the availability of a large number of unclassified OTUs at the genus level in 16S databases, as well as the ability of amplicon-based sequencing to capture low-abundance taxa. In addition, that a number of observed genera were not identified using the marker gene-mapping-or the alignment-based approach indicates that reference databases will continue to benefit from the sequencing of additional novel bacterial. The ability to detect low-abundance bacterial genera and species being impaired at shallow sequencing depths when using the marker gene-mapping-based approach suggests that the combination of shallow sequencing and marker gene-mapping-base annotation may not be suitable for studies seeking to characterize low-abundance bacterial taxa, including those that are transient members of the gut microbiome. Notably, the number of bacterial genera and species detected were greater when using the alignment method compared to both the marker gene-mapping-based approach and 16S sequencing. This may be in part due to differing numbers of reference genomes in the interrogated databases, or alignment of reads to closely related taxa. While only one 16S database and variable region were tested, results open the opportunity to evaluate the performance of additional 16S databases, as well as other 16S variable regions in comparison with SMS results. α-diversity and relative abundance results from the selected taxa showed that the alignment-based method may have a greater concordance in signal across sequencing depths and indicate that this approach may be more suitable when using SMS with the caveat of having reads aligned to closely related bacteria, which could result in false positives ( Table 2). In addition, these data indicate that, in some cases, 16S sequencing may be more suitable depending on the organism of interest. These data demonstrate that the selection of both sequencing depth and annotation method for taxonomic classification should be carefully evaluated and, ideally, include reference materials to assess performance. Results from reference materials in the present study demonstrated that further filtering of data may be needed to decrease the number of false positives with the caveat of losing taxa present in extremely low relative abundances. This was the case for the alignment-based approach where the selected filtering percentage (i.e., 0.01%) reduced the number of false positives. Results suggest that further filtering may need to be tested accordingly. Table 2. Summary of shotgun metagenomic sequencing results at various sequencing depths and annotation with marker gene mapping (MetaPhlAn2) and alignment (BURST) at genus and species levels. Genus level results were compared to 16S information. Table also shows summary of pathway information annotated using UniProt and MetaCyc, as well as virome information at various sequencing depths.

Marker Gene-Mapping
• α-diversity: Number of identified genera is slightly lower compared to alignment and 16S methods. May require > 0.5 Gb to capture diversity similar to 5 Gb. • β-diversity: Grouping of the data is subject-and not sequencing-depth-based. • Taxonomy/Abundance: Detection may depend on relative abundance (i.e., low abundant bacteria may not be detected at lower sequencing depths).

Alignment
• α-diversity: Number of identified genera is higher compared to marker gene-mapping and 16S methods. May require > 0.25 Gb to capture diversity similar to 5 Gb. • β-diversity: Grouping of the data is subject-and not sequencing-depth-based.

•
Taxonomy/Abundance: Detection may not depend on relative abundance (i.e., bacteria in low abundances may be detected at lower sequencing depths).

16S
• α-diversity: Number of identified genera is slightly higher than marker gene-mapping, but lower than alignment approach with additional filtering of taxa present at relative abundances < 0.01%. • β-diversity: Data group more closely to alignment data. • Taxonomy/Abundance: Genus level resolution achieved.

Marker gene-Mapping
• α-diversity: Number of identified species is slightly lower compared to alignment method. May require > 0.5 Gb to capture diversity similar to 5 Gb. • β-diversity: Grouping of the data is subject-and not sequencing-depth-based. • Taxonomy/Abundance: Detection may depend on relative abundance (i.e., low abundant bacteria may not be detected at lower sequencing depths).

Alignment
• α-diversity: Number of identified species is slightly higher compared to marker gene-mapping method. May require > 0.5 Gb to capture diversity similar to 5 Gb. • β-diversity: Grouping of the data is subject-and not sequencing-depth-based. • Taxonomy/Abundance: Detection may not depend on relative abundance (i.e., bacteria in low abundances are detected at lower sequencing depths). Alignment to closely related taxa may result in false positives.

16S
• One variable region (i.e., V4) did not discriminate species level. SMS also exerts an effect on the lower-abundance, less-characterized aspects of the microbiome. By using viruses as an example, this study demonstrates that the effect of the sequencing depth is more profound when studying the less-studied members of the human microbiome. While the annotation method may also affect virus detection and identification, testing different tools was outside the scope of the present study and has previously been examined in comparison with VirMap [28]. Although not all viruses were identified in similar relative abundances in all subjects, as viromes tend to be personalized, results here provide proof-of-concept data covering various subject-to-subject and virus-to-virus scenarios at various sequencing depths. While we were limited to 5 Gb, it seemed that > 5 Gb may be best suited to capture DNA viruses from SMS data when no a priori isolation/concentration steps for viruses have been applied (Table 2). Results may also translate to other low abundant categories of the microbiome, including potential biomarkers, virulence factors, and antibiotic-resistant genes.
While taxonomic assignments provide key information in microbiome studies, pathway information has the ability to provide insights into the potential role(s) of specific microbial members in maintaining health, promoting disease, and/or be associated/correlated with specific processes. Understanding pathway profiles also represents an opportunity to develop metabolite-based therapeutics, also known as "postbiotics" [34]. Various well-known examples of postbiotics include short-chain fatty acids (SCFAs) [35], taurine [36], and flavonoids [37]. SMS provides insights into the identification of genes associated with the metabolism of some of these specific molecules, but sequencing depth also has the potential to influence information recovery and data interpretation. SSMS may be suitable for a global pathway profile characterization, but it may not be suitable for the characterization of rare pathways, or those that are present at low relative abundances. For this reason, the expected prevalence and abundance of each pathway must be considered to determine the specific sequencing depth needed. Nevertheless, a sequencing depth of approximately 1 Gb in human stool samples appears to stably capture a large proportion of pathways. Future studies need to focus on determining the performance of various pathway annotation tools at various sequencing depths and sample types using appropriate reference materials, as needed (Table 2).

Conclusions
In this study, we evaluated the combined effects of depth of sequencing and annotation method (i.e., marker gene mapping-vs. alignment-based approaches) on information recovery for SSMS. Sequencing depths > 0.25 Gb and > 0.5 Gb capture α-diversity and genus-level taxonomic abundances similar to those recovered at sequencing depth of 5 Gb when using marker gene-mapping and alignment methods, respectively. Sequencing depths > 0.5 Gb captured bacterial species information at levels analogous to those detected with a sequencing depth of 5 Gb. Pathway recovery at sequencing depths > 0.75 Gb was comparable to that achieved at 5 Gb. Virus profiles were individualized, and signal recovery was impaired at sequencing depths < 5 Gb. Results suggest that sequencing depth and annotation method need to be carefully evaluated depending on the aims of the study.