Impact of Preservation Method and 16S rRNA Hypervariable Region on Gut Microbiota Profiling

Large-scale human microbiota studies require specimens collected from multiple sites and/or time points to maximize detection of the small effects in microbe-host interactions. However, batch biases caused by experimental protocols, such as sample collection, massively parallel sequencing, and bioinformatics analyses, remain critical and should be minimized. This work evaluated the effects of preservation solutions and bacterial 16S rRNA gene primer pairs in revealing human gut microbiota composition. Since notable changes in detecting bacterial composition and abundance were observed among choice of preservatives and primer pairs, a consistent methodology is essential in minimizing their effects to facilitate comparisons between data sets.

associated microbes reside, indicate that the disruption of microorganisms and their abundances could be linked to a number of diseases, such as carcinogenesis (3), cardiovascular disease (4), diabetes (5), hypertension (6), inflammatory bowel disease (7), and obesity (8). The study of gut microbiome has been empowered in a cultureindependent manner through high-throughput sequencing and bioinformatics analyses (9)(10)(11). Because composition of gut microbiota differs greatly between individuals, microbiome studies have expanded into larger populations at multiple sites and/or time points to maximize detection of the small effects of microbe-host interactions. However, biases caused by methods of sample collection and choice of 16S rRNA gene hypervariable region remain critical and should be minimized.
Immediately freezing stool samples at Ϫ20°C or below is considered the "gold standard" to prevent shifting of microbial community composition (12). However, this approach is not feasible for sample collection in remote areas without reliable cold chain transport, or in studies where subjects are requested to send self-collected household samples to laboratories at ambient temperature. The effects of several storage conditions in revealing microbiota composition have been compared to immediate freezing at Ϫ80°C. For example, stool samples stored at room temperature for at most 1 day, or at 4°C and Ϫ20°C for up to 2 weeks, had little effect on microbiota composition shift (13)(14)(15). When samples are exposed to room temperature or higher other than ultralow freezing, however, use of preservation methods becomes essential and critical. Ethanol in concentrations of 95% or higher preserved bacterial DNA for long-term storage well (15), but another study has reported low DNA yield as a potential drawback of using ethanol (16). In addition, self-collected household samples stabilized with ethanol (70% and above) might require special shipping, depending upon local transportation regulations. RNAlater, commonly used as a general preservative for RNA samples, has been shown to have decreased DNA purity and to perform relatively poorly in maintaining microbiota composition derived by 16S rRNA gene amplicon sequencing (16,17). Buffers containing EDTA can inhibit the growth of certain bacteria (18), but they significantly influenced gut microbiota profiles, such as an increased abundance of Bacteroides and Proteobacteria and a reduction of Firmicutes and Actinobacteria (13). OMNIgene.Gut, a commercially available approach to preserve stool samples at ambient temperature, has recently been reported to generate community profiles that diverge the least from Ϫ80°C standards (13,19). However, its effectiveness relative to other approaches needs to be independently assessed over longer time scales. Moreover, these prior studies have generally not evaluated the bactericidal ability of preservatives in rendering microbiota samples noninfectious.
Since microbiome research is quickly shifting toward determination of small effect size variations relevant to host health, understanding the effects of preservative solutions is crucial to minimize bias. As part of the protocol development for a multiple-site, large-scale longitudinal microbiota study, we evaluated the effects of six preservative solutions (Norgen, OMNI, RNAlater, CURNA, HEMA, and Shield) in retaining gut microbiota composition compared to immediate freezing at Ϫ80°C. Stool samples from five healthy individuals were collected and subjected to storage in different preservation solutions for 7 days, a time period typically sufficient for delivery to the laboratory by post. Three hypervariable regions of the bacterial 16S rRNA gene (V1-V2, V3-V4, and V4) were amplicon sequenced. The bactericidal ability of these preservatives was assessed by standard aerobic and anaerobic cultures.

RESULTS
Study subjects. Five healthy individuals provided fresh stool samples ( Fig. 1; see also Table S1 in the supplemental material). For each subject, one subset of sample was immediately frozen at Ϫ80°C; eight aliquots were stored in different solutions, including six preservatives (Norgen, OMNI, RNAlater, CURNA, HEMA, and Shield), phosphatebuffered saline (PBS) solution, and Cary-Blair (CB) transport medium (Table 1). Following 7 days of preservation at room temperature, all aliquots were stored at Ϫ80°C without a freeze-thaw cycle until DNA extraction. Total DNA, including Ϫ80°C "gold standards," was extracted with a bead-beating method and amplicon sequenced on an Illumina MiSeq platform by targeting three hypervariable regions of the bacterial 16S rRNA gene (V1-V2, V3-V4, and V4) ( Table 2). No obvious difference in quality of extracted DNA for PCR amplification was observed between subjects or preservatives, as visualized on the agarose gel with similar band brightness (data not shown).
Variation in gut microbial communities without preservative at room temperature. We first used sequences representing the V3-V4 region of the 16S rRNA gene to assess gut microbiota profile since this region provides the longest read length ( Table 2). A total of 121,532 high-quality sequence reads for 45 samples (mean Ϯ SD of 2,700 Ϯ 667) was obtained, clustering into 308 amplicon sequence variants (ASVs). These ASVs represented 105 bacterial taxa at genera or higher levels, while 64 were detected at Ն1% relative abundance in at least one sample. The gut microbial community assigned at order level is displayed in Fig. 2A. As expected, the nonpreserved FIG 1 Workflow of the study design. Nine stool sample aliquots from each of the five donors were subjected to different preservatives and incubated at room temperature for 7 days to compare shifts in gut microbial community composition. stool samples stored in PBS (8_PBS) and Cary-Blair transport medium (9_CB) at room temperature for 7 days had significant changes in microbiota composition compared to Ϫ80°C standards and preserved stocks (samples 1 to 7), as discriminated by a principalcoordinate analysis (PCoA) using either weighted or unweighted UniFrac distances ( Fig. 2B and C). Using the Adonis function in R's package Vegan, however, a permutational multivariate analysis of variance (PERMANOVA) based on weighted UniFrac distances (df ϭ 2, R 2 ϭ 0.504, pseudo-F ϭ 21.346, P Ͻ 0.001) found a higher variability in microbial community composition than unweighted UniFrac distances (df ϭ 2, R 2 ϭ 0.093, pseudo-F ϭ 2.160, P Ͻ 0.01) in separating nonpreserved samples from preserved ones, indicating that the sample clustering was driven more by the proportion of microbial community members than the presence/absence of bacterial taxa. For example, a significant increase of the relative abundance of Enterobacteriales (Fig. 2D) and Fusobacteriales (Fig. 2E) and a reduction of Clostridiales (Fig. 2F) and Betaproteobacteriales ( Fig. 2G) were observed in nonpreserved gut samples (Table S2). Variation in gut microbial communities stored in preservatives. Next, we refined the ASV table by excluding samples in PBS solution and Cary-Blair transport medium to measure the variation of six preservatives in retaining gut microbiota composition with reference to immediate freezing at Ϫ80°C as standards. Overall, the microbial community of the surveyed stool samples was mainly dominated by Faecalibacterium (mean relative abundance of 23.4%) and Bacteroides (22.2%), followed by Roseburia (7.3%) ( Table S3). The alpha diversity analyses at either ASV or genus level did not differ significantly between Ϫ80°C standards and preserved samples ( Fig. 3A and B), except that aliquots stored in Shield appeared relatively lower in Shannon diversity (mean of 3.0 and 2.2, respectively) and Simpson evenness (mean of 0.9 and 0.8, respectively). We found significant changes of the alpha diversities of gut microbial community across subjects (Kruskal-Wallis test, P Ͻ 0.029) (Fig. S1).
A PERMANOVA using weighted UniFrac distances indicated that both subject groups (df ϭ 4, R 2 ϭ 0.632, pseudo-F ϭ 26.134, P Ͻ 0.001) and preservation conditions (df ϭ 6, R 2 ϭ 0.223, pseudo-F ϭ 6.157, P Ͻ 0.001) contributed to the differences in microbiota composition, although the majority of variability was from differences between subjects rather than the lower impact of storage conditions (63.2% versus 22.3%). The variation was statistically attributed to subjects (df ϭ 4, R 2 ϭ 0.779, pseudo-F ϭ 28.034, P Ͻ 0.001) but not preservatives (df ϭ 6, R 2 ϭ 0.055, pseudo-F ϭ 1.309, P ϭ 0.118) based on unweighted UniFrac distances. In line with these observations, the principal-coordinate analysis supported that microbial variation was mainly associated with differences between individuals ( Fig. 3C and D). In order to examine the small effects of variation by preservation methods, we subsequently calculated the pairwise distances between Ϫ80°C standards and preserved samples collected from the same subjects using weighted UniFrac distances, since unweighted UniFrac distances did not reveal statistical significance (P ϭ 0.997) (Fig. 3C). As shown in Fig. 3D (weighted), samples stored in various preservatives exhibited diversified microbial community shifts relative to Ϫ80°C standards (intrasubject distance to standard of 0.06 to 0.26) (P ϭ 0.027), suggesting that all preservation methods contributed impacts on the variation in microbiota profiling. We found, however, that two preservatives (Norgen and OMNI) showed the least shift in community composition (mean distance to the standards of 0.09 and 0.09, respectively) while RNAlater had the highest impact (mean samples stored in Shield (P ϭ 0.0239) and a decrease of the relative abundance of Lachnospiraceae (Eubacterium eligens group) in OMNI (P ϭ 0.0461) relative to Ϫ80°C standards were observed. Impact of targeting different regions of the 16S rRNA gene in profiling gut microbiota. Previous work has suggested that targeting different variable regions of the 16S rRNA gene may generate different microbial community profiles (20). In this study, we compared the different performances of primer pairs targeting three hypervariable regions of bacterial 16S rRNA gene (V1-V2, V3-V4, and V4) in profiling gut microbiota composition. Overall, 254,363 and 420,622 high-quality sequence reads (mean Ϯ SD of 7,267 Ϯ 1,139 and 12,017 Ϯ 3,799) were obtained with V1-V2 and V4 primer pairs, respectively. When analyses on the divergences of preservation methods in revealing gut microbial community were repeated, the sequences generated by V1-V2 and V4 primer pairs ( Fig. 4; Fig. S2 and S3) showed similar results as that of V3-V4 reads. The surveyed stool samples mainly clustered by subjects. Although variations in microbiota composition were observed for samples in all preservatives tested, the effects were relatively small (PERMANOVAs using weighted UniFrac distances: V1-V2, R 2 ϭ 0.163, P Ͻ 0.001; V4, R 2 ϭ 0.167, P Ͻ 0.001) compared to the high heterogeneity between individuals (V1-V2, R 2 ϭ 0.709, P Ͻ 0.001; V4, R 2 ϭ 0.709, P Ͻ 0.001). Samples in Norgen preservative consistently showed the least shift in community composition relative to Ϫ80°C standards. Meantime, the relative abundance of Faecalibacterium in samples stored in Shield was significantly higher, as revealed by all three 16S rRNA gene region primer pairs (Table S3). However, the choice of primers targeting different hypervariable regions of the bacterial 16S rRNA gene may specifically impact microbial community profile and bacterial taxon assessment. When all three region reads were combined for principalcoordinate analyses and PERMANOVAs using fragment insertion SATé-enabled phylogenetic placement (SEPP) technique, we found strong bias of gene regions in sample clustering based on either unweighted (df ϭ 2, R 2 ϭ 0.218, pseudo-F ϭ 42.239, P Ͻ 0.001) or weighted (df ϭ 2, R 2 ϭ 0.197, pseudo-F ϭ 51.940, P Ͻ 0.001) UniFrac distances, although the majority of variability was consistently from differences between subjects (unweighted, df ϭ 4, R 2 ϭ 0.527, pseudo-F ϭ 51.073, P Ͻ 0.001; weighted, df ϭ 4, R 2 ϭ 0.507, pseudo-F ϭ 66.863, P Ͻ 0.001) (Fig. 5A and B). The impact of preservation methods on microbial community shifts was relatively small (unweighted, df ϭ 6, R 2 ϭ 0.018, pseudo-F ϭ 1.156, P ϭ 0.198; weighted, df ϭ 6, R 2 ϭ 0.122, pseudo-F ϭ 10.707, P Ͻ 0.001). We also found that the V4 prime pair showed a higher Shannon  Table S4). diversity for gut microbiota profile than did other primer sets (Fig. 5C), in line with the higher levels of alpha diversity with rarefied V4 region reads (at genus level and above) (Fig. S4). The relative abundance of 11 bacterial orders was significantly variable among reads generated by different primers (Fig. 5D; Table S4). For example, the commonly used V1-V2 primer pair failed to detect the majority of Bifidobacteriales (mean relative abundance of 0.03% versus 6.40%, P Ͻ 0.001) but revealed more Enterobacteriales (5.69% versus 2.57%, P Ͻ 0.001) and Erysipelotrichales (3.2% versus 2.0%, P Ͻ 0.001) compared to another two primer sets. The V3-V4 primers showed the highest abundance of Clostridiales (56.69% versus 44.12%, P Ͻ 0.001), which likely led to the apparent decrease of Bacteroidales (26.91% versus 34.97%, P Ͻ 0.001), Betaproteobacteriales (2.45% versus 3.58%, P Ͻ 0.001), Coriobacteriales (0.26% versus 0.71%, P Ͻ 0.001), and Pasteurellales (0.13% versus 0.31%, P Ͻ 0.001). Mollicutes RF39 was amplified only by the V4 primer set, although the relative abundance of this bacterial group was extremely low in the surveyed stool samples.
Sequence variation of primer-binding sites. In order to characterize the occurrences of sequence variations at the primer-binding sites used in this study, we refined an NCBI bacterial 16S rRNA reference database (a total of 18,773 complete or nearcomplete sequences) using a threshold of 99% similarity and extracted sequence sets matching each primer-binding site. The most common sequence variations are shown in Fig. 6, with the corresponding dominant bacterial group(s) in which the sequence variations are observed. It is noted that the V1-V2 forward primer-binding site sequences in most Bifidobacteriales and Chlamydiales 16S rRNA genes were not precisely accommodated by the commonly used degenerate 27f primers, differing at two and three positions from the best-matching 27f-YM primer, respectively. Similarly, the nondegenerate form of the V1-V2 reverse primer used in this study (338R) may have low sensitivity in detecting Planctomycetales, for which three mispairings were identified. The V1-V2 reverse primer (338R) and V3-V4 forward primer (341F) showed mismatches with Mollicutes RF39 16S rRNA genes in two positions, which could confer a low efficiency in detecting this group of bacteria (Fig. 5D). Optimization of primer pools to minimize the mispairing at primer-binding sites is proposed (Fig. 6B).
Bacterial culture of stool samples in different preservatives. To test whether the surveyed preservatives are able to rend microbes noninfectious, stool suspensions preserved in each solution over time were cultured for bacteria under both aerobic and anaerobic conditions (Fig. S5). We found that Norgen, OMNI, HEMA, and Shield preservatives immediately and completely prevented the growth of gut bacteria, allowing the resulting samples to be handled safely. In contrast, RNAlater and CURNA, two preservation solutions designed for RNA samples, as well as the PBS and Cary-Blair transport medium, showed no effect in inactivating bacteria.

DISCUSSION
Microbiota studies have been quickly switched to large-scale sample collection to characterize the subtle differences in microbial communities in relation to human health and diseases. Although immediate freezing at Ϫ20°C or below has been considered the "gold standard" for microbiota preservation, this approach is not feasible for many studies, especially field sampling from remote areas or self-collected samples requiring shipment to the laboratory without cold chain transport. In such cases, a proper preservation method to minimize microbiota composition shift to reduce bias and to inactivate infectious agents for samples to be handled safely becomes critical. Hence, we evaluated the effects of multiple preservation solutions and 16S rRNA gene hypervariable regions in revealing gut microbiota composition. Our data showed that the profile of gut microbiota can change dramatically over a course of 7 days at room temperature in the absence of protection of preservation. The use of a proper preservative without cold chain, in contrast, provides satisfactory prevention of community shifts. The effects of these preservatives were small compared to intersubject variations; however, notable changes in microbiota composition were observed. Nevertheless, we found that two preservatives, Norgen and OMNI, showed the smallest microbial community shifts, and both efficiently inhibited the growth of aerobic and anaerobic bacteria, suggesting they may serve as useful solutions for stool sample transportation and storage. Recent studies have reported that OMNIgene.GUT (here called OMNI) is able to stabilize DNA at room temperature for up to 14 days, with little difference in microbiota composition by 16S rRNA gene sequencing (13,15,19,21). RNAlater, a commonly used preservative for RNA samples, on the other hand, exhibited larger community shift in bacterial relative abundance at room temperature for 7 days, consistent with previous reports that the microbial community in RNAlater may lose stability at longer time frames when maintained at ambient temperatures (12,15,16,22). However, under cold conditions (e.g., 4°C or Ϫ20°C) or for a short period of room-temperature storage (e.g., within same or another day), it performed equally well as other storage conditions, such as 95% ethanol, FTA cards, and OMNIgene.GUT (21,23,24). Interestingly, RNAlater and CURNA (with similar components as RNAlater) do not prevent bacterial activity, highlighting themselves as potential alternatives for microbial sample collection when live bacterial culture is necessary. According to the manufacturer's instructions, HEMA was initially developed for the stabilization of buffy coat samples and there was no previous report evaluating its use as preservative for samples intended for microbiota studies. Shield performed well in microbiota preservation and disinfection but consistently showed increased relative abundance of Faecalibacterium regardless of the choice of 16S rRNA gene primer pairs. We recommend the use of Norgen or OMNI for microbial sample preservation and caution against RNAlater for long-term storage at ambient temperatures.
Choosing primers to target certain regions of the 16S rRNA gene for bacterial taxonomic profiling is another difficult challenge, and the comparability of primers with microbes from various communities requires attention (25,26). The V1 to V4 region, for example, is more divergent than V4 to V7 and can provide a higher resolution in bacterial clustering for stool samples (20,27). In this work, we compared the performance of three commonly used primer pairs targeting the 16S rRNA gene V1-V2, V3-V4, and V4 regions in profiling gut microbial community and found a higher alpha diversity and richness revealed by V4 primers. Our results also indicate an immense impact of primers with regard to amplification of bacterial taxa. Close examination of primerbinding sites revealed several sequence variations within cohesive phylogenetic groups (e.g., Bifidobacteriales) that are not accommodated by the commonly used 27f primer sequence. In such cases, the use of pooled primers to minimize mismatches with this bacterial group will provide an effective approach for revealing the true bacterial diversity (26,28,29). Knowledge of how primer pairs differentially amplify bacterial taxa is important for study design and relating results to prior studies. Alternatively, the quickly developing next-generation sequencing platforms, such as PacBio and Nanopore, may facilitate the characterization of full-length bacterial 16S rRNA genes and allow global comparisons of microbiome studies across sample types and gene regions (30)(31)(32).
This study has its strengths and limitations. We expand the current understanding of the performance of several DNA and RNA preservation solutions, including CURNA and HEMA, to our knowledge being evaluated for the first time, in retaining gut microbiota composition. The data provide a useful reference for knowing the biases introduced by storage conditions. We have expanded the knowledge of bactericidal ability of the surveyed preservatives that may provide alternatives in study design considering the potential of infectious agents and bacterial culture. We have also identified the significant difference of 16S rRNA gene primer pairs in charactering bacterial taxa and emphasized the importance of a consistent methodology for large epidemiological studies. Nevertheless, one of the shortcomings of this study is the small sample size and the conserved sample source (5 subjects who were all male, young, and healthy), which may underestimate the impact of preservation solutions and primer pairs in charactering the diversity of gut microbiota. Further methodological comparison is needed to expand to other settings, including patient groups and population-based epidemiological studies. Only sequencing-based evaluation was per-formed in this work, while the concordance of preservation solutions in other omics studies, such as proteomics and metabolomics, may be largely variable. For example, RNAlater was reported to be not feasible for metabolomics measures probably due to its high sodium sulfate content, making the preserved samples incompatible with metabolomics platforms (22,33). Last, more robust qualitative and quantitative testing under variable conditions, such as extreme temperatures, increased freeze-thaw cycles, or longer storage periods, may provide a more comprehensive evaluation of the effectiveness and stabilization of storage methods at preserving microbial DNA quality and yield.
Conclusions. It is important for large-scale microbiome studies to accurately and consistently reveal microbial communities while minimizing external effects. Choices of preservation solution for transportation of stool samples at ambient temperature and of gene region for 16S rRNA gene amplicon sequencing are critical determinants impacting the observed microbial community profiling. We recommend the use of a proper preservative, such as Norgen or OMNI, that can minimize microbial community shift and efficiently inactivate bacterial growth, allowing the resulting samples to be handled and shipped safely and stably. Amplicon sequencing-based human microbiota studies are highly dependent on choice of primer pairs targeting different hypervariable regions of the bacterial 16S rRNA gene. Nevertheless, a consistent methodology is essential in minimizing batch biases to facilitate comparisons between data sets.

MATERIALS AND METHODS
Ethics approval. This study was approved by the Joint Chinese University of Hong Kong-New Territories East Cluster Clinical Research Ethics Committee. All subjects recruited in this study were older than 18 years of age, and samples were anonymized without individual identifying information. Written informed consent was obtained from each participant.
Stool specimen collection, processing, and DNA extraction. Five healthy male subjects, aged between 23 and 28 years old, were recruited to provide fresh stool samples using a stool self-collection kit (see Table S1 in the supplemental material). Stools from these five donors were collected on different days and were transported to the laboratory in an icebox cooler at 4°C within 2 h without a freeze-thaw cycle. Upon receipt, approximately 5 g of sample was immediately mixed with 200 l of sterile phosphate-buffered saline (PBS) and thoroughly homogenized for 1 min using a sterile cotton stick. One gram of the homogenized stool, serving as a standard, was transferred into a 2-ml sterile collection tube and immediately stored at Ϫ80°C until further DNA extraction. In parallel, eight stool aliquots, each weighing approximately 0.2 g, were transferred into 2-ml collection tubes preloaded with 1 ml of preservative (Norgen, OMNI, RNAlater, CURNA, HEMA, and Shield), PBS, or Cary-Blair transport medium, which is intended for preserving the viability of enteric bacteria (Table 1). Tubes were pulse vortexed for 5 s and placed on a laboratory bench at room temperature for 7 days. Then, stool aliquots were stored at Ϫ80°C without a freeze-thaw cycle for 28 to 73 days until DNA extraction (Fig. 1). The procedures of stool homogenization and aliquoting were performed in a biosafety cabinet and finished within 30 min.
DNA extraction was performed in accordance with the Earth Microbiome Project protocol (34), with the bead-beating QIAamp PowerFecal DNA kit (Qiagen, USA). In detail, 100 l of stool suspension from aliquots stored at Ϫ80°C was added to the PowerBead tubes for total DNA extraction according to the manufacturer's protocol. An MP FastPrep-24 5G homogenizer, at a speed of 6.0 m/s for 45 s, was used for mechanical lysis. Approximately 0.2 g of stool standards was suspended in 1 ml of sterile PBS prior to DNA extraction. The DNA samples were eluted in 100 l elution buffer (pH 8.0) and stored at Ϫ20°C. DNA concentrations were measured using a NanoDrop 2000 (Thermo Scientific, USA).
16S rRNA gene PCR amplification and sequencing. We used three PCR primer pairs targeting the 16S rRNA gene V1-V2, V3-V4, and V4 regions, respectively, to profile gut microbiota compositions ( Table 2) (35). A pair of dual 12-bp barcodes was indexed to each amplicon set through the forward and reverse primers modified from the Earth Microbiome Project protocol (34). Successful amplicons were pooled at approximately equal molar DNA concentrations, purified using the QIAquick gel extraction kit (Qiagen, USA), and sequenced on an Illumina MiSeq (Illumina, USA) at the Weill Cornell Medicine Genomics Resources Core Facility, New York, NY, USA, using paired-end 300-bp reads.
Bioinformatics and biostatistics analysis. Following demultiplexing, the QIIME2 (2018.8) package (https://qiime2.org) (36) was applied to assign Illumina short reads into amplicon sequence variant (ASV) and bacterial taxa tables. The detailed pipelines, including quality control, paired-end short reads merging, and dada2 denoising and clustering, can be found in the supplemental pipeline (Text S1). In order to combine reads from different 16S rRNA hypervariable regions for principal coordinates analysis, the fragment insertion SATé-enabled phylogenetic placement (SEPP) technique with default parameters was used (37). The SILVA 132 99% 16S rRNA gene reference database (https://www.arb-silva.de/ download/archive/qiime) was used to assign bacterial taxonomic classification (38). Singleton reads were removed for statistical analysis. Compositions of microbiota communities were summarized by proportion at different taxonomy levels, including genus, family, order, class, and phylum ranks.
In order to describe the diversity of observed bacterial taxa taking into consideration richness and abundance, Shannon and Simpson diversity indexes (39) were calculated. We assessed the alpha diversity using species richness and Shannon diversity at various rarefaction depths ranging from 1 to 6,001 reads (with steps of 250 reads) and found no significant difference of alpha diversity using rarefied data sets at depth of 1,000 reads per sample or higher (Fig. S4). In order to retain all samples for diversity analysis, we set 1,500 reads per sample as rarefaction depth to normalize the data for differences in sequence count. We used a UniFrac algorithm in the GUniFrac R package (40), with unweight or weight (alpha ϭ 1) on abundant lineages to calculate pairwise Kantorovich-Rubinstein (KR) distances between samples. Differences in community composition were assessed using permutational multivariate analysis of variance (PERMANOVA) in the Vegan R package. Principalcoordinate analysis was performed to visualize associations between community composition and experimental factors. Comparisons of the relative abundances of characteristic ASVs between defined groups were performed using nonparametric Mann-Whitney Wilcoxon rank sum test (MWU), Wilcoxon signed rank test (WSR), Kruskal-Wallis test (KW), or Tukey's honest significant difference (Tukey HSD) post hoc test. Statistical analyses and plotting were performed in R (3.4.0). A two-sided P value of Ͻ0.05 was considered statistically significant.
Bacterial aerobic and anaerobic culture. Approximately 0.2 g of homogenized fresh stool samples was transferred into 2-ml collection tubes preloaded with 1 ml of preservative (Norgen, OMNI, RNAlater, CURNA, HEMA, and Shield), PBS, or Cary-Blair transport medium. For each sample, 100-l replicate suspensions were incubated at room temperature for 0 min, 10 min, 1 h, 5 h, 24 h, 48 h, and 96 h, respectively. Approximately 100 l of 1:10-diluted stool suspensions in 0.85% saline solution was spread on a blood agar plate and a K1 agar plate for aerobic and anaerobic culture, respectively. Plates were incubated at 37°C for 48 h in a duplicate setting, and the number of CFU was counted immediately afterwards.
Data availability. All sequence data generated from this study were deposited in the Sequence Read Archive with accession number PRJNA470603.