The impact of storage conditions on human stool 16S rRNA microbiome composition and diversity

Ethics statement

This study was approved by the University of Glasgow College of Medical Veterinary and Life Sciences Ethics Committee (project code 200160068), the Vector Control Division, Ministry of Health Uganda, Research Ethics Committee (reference: VCDREC/062) and the Uganda National Council for Science and Technology (UNCST-HS 2193). Informed signed or fingerprinted parental or guardian consent, and signed or fingerprinted assent from the study children was obtained prior to participation.

Sample collection

One stool specimen was collected from three children, aged 12–14, selected at random from Bugoto Lake View Primary School, Mayuge District, Uganda in March 2017. The sample from Child A was collected on day 1, and those from Child B and Child C on day 2. Outer surface, central inner and mixed stool samples (~300 mg each) were taken from each specimen and stored separately in cryovials as raw stool (considered the standard), dispersed in absolute ethanol (approx. stool:ethanol ratio = 1:6) and dispersed in RNAlater (approx. stool:RNAlater ratio = 1:6), then frozen immediately on dry ice. Additionally coarsely homogenised stool from each donor were frozen on dry ice at 1, 2, 4, 8, 16 and 32 h post collection for each storage preservation method (Table S1). Time zero was taken as the time at which all the stool samples, taken from an individual stool specimen, had been processed into all the collection tubes for the relevant conditions to be tested, which was approximately 30 min after defecation. Prior to freezing on dry ice, stool was kept in cooler, shaded, well ventilated, indoor spaces as much as realistically possible. Within 48 h of freezing on dry ice, samples were transferred into a −20 °C freezer and later transported to the University of Glasgow on dry ice for further processing and analysis. Samples underwent one freeze-thaw cycle (<30 min) during weighing and, to the best of our knowledge, they remained frozen at −20 °C from collection until DNA was extracted approximately 6 months later. Cryovials used in the field containing only ethanol or RNAlater, without stool, were used as negative controls. Samples from two of the children were also stored using OMNIgene.GUT kits (DNA Genotek (Doukhanine et al., 2014)) as per manufacturer’s instructions, and remained at ambient temperature until DNA was extracted approximately six (four in Uganda and two in the UK) months later (~three times the recommended 60 day stability recommendation for the kit; http://www.dnagenotek.com/us/products/collection-microbiome/omnigene-gut/OMR-200.html). Details of samples and sample codes are shown in Table S1.

Extraction of DNA from stool

The MPbio FastDNA™ SPIN Kit for Soil (MP Biomedical, Irvine, CA, USA), was used to extract nucleic acids from ~200 mg of stool with minor modifications to the method described by Alcon-Giner et al. (2017), as follows. An attempt was made to exclude large pieces of undigested vegetable matter from stool during the weighing process. Samples were homogenised using a TissueLyser II (Qiagen, Hilden, Germany) at a speed setting of 25, and a 2 min centrifugation was used after addition of binding matrix. DNA concentration was quantified using a NanoDrop 1000 fluorimeter (Thermo Fisher Scientific, Waltham, MA, USA).

16S library preparation

A modification of the Illumina 16S metagenomic sequencing library preparation protocol (https://support.illumina.com/documents/documentation/chemistry_documentation/16s/16s-metagenomic-library-prep-guide-15044223-b.pdf) was used to prepare the DNA library. PCR was used to amplify the V1–V2 regions of the 16S rRNA gene, chosen because they were better at detecting bacterial species of interest from stool for future studies (eg. Bifidobacterium (Alcon-Giner et al., 2017)). The primers used were: 16SV1 forward primer (5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGAGMGTTYGATYMTGGCTCAG-3′) and 16SV2 reverse primer (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGCTGCCTCCCGTAGGAGT-3′).

Each reaction was performed in a final volume of 25 µL consisting of 1× KAPA HiFi HotStart ReadyMix (KAPA Biosystems, Wilmington, MA, USA), 0.5 µM of each primer and 12.5 ng of sample DNA. Thermocycler conditions were used as follows: 95 °C for 5 min, followed by 26 cycles of 95 °C for 30 s and 60 °C for 1 min. Samples were then held at 10 °C in the PCR machine, before being stored at 4 °C. H₂O sample controls were included as negative controls during the first round of PCR to monitor non-specific amplification.

Each PCR product was purified by mixing with a 0.90 × PCR product volume of High Prep PCR beads (MAGBIO, Gaithersburg, MD, USA). After a 10 min incubation, sample tubes were placed on a magnetic stand and left until the supernatant became clear. The supernatant was then removed and the beads were washed twice with freshly prepared 80% ethanol, and then left to dry for 15 min to allow residual ethanol to evaporate. The sample tubes were removed from the stand and the beads were then resuspended in 20 µL Tris buffer pH 8.5, and incubated for 2 min before being placed back on the magnetic stand. Once clear, the supernatant was transferred to a fresh tube and the DNA concentration quantified using the Quant-iT PicoGreen dsDNA Assay (https://assets.thermofisher.com/TFS-Assets/LSG/manuals/mp07581.pdf) (Thermo Fisher Scientific, Waltham, MA, USA).

A second PCR step was then used to barcode each sample. PCR reactions were performed in a final volume of 50 µL consisting of 1× KAPA HiFi HotStart ReadyMix, 5 µL of each of two Nextera XT Index Kit Set A (Illumina, San Diego, CA, USA) indices, with each sample having a unique combination, and 10 ng of post-PCR1 sample DNA. Thermocycler conditions used were as follows: 95 °C for 3 min; followed by eight cycles of 95, 55 and 72 °C for 30 s each; with a final step of 72 °C for 5 min. Samples were then held at 10 °C in the PCR machine, before being stored at 4 °C.

Samples were cleaned with High Prep PCR beads as described above and then combined to form an equimolar sample library. The Wizard SV Gel and PCR Clean-Up System Kit (Promega, Madison, WI, USA) was used to purify the DNA library prior to sequencing, as per manufacturer’s instructions (https://www.promega.co.uk/-/media/files/resources/protocols/technical-bulletins/101/wizard-sv-gel-and-pcr-clean-up-system-protocol.pdf?la=en) using a band size of ~435 bp. DNA concentration was then measured using a Bioanalyser 2100 (Agilent, Santa Clara, CA, USA).

Sample sequencing and analysis

Samples were sequenced using the Illumina MiSeq platform (Glasgow Polyomics, Glasgow, UK) with 2 × 300 bp paired-end read lengths with up to 100,000 reads per sample (MiSeq V3 600 cycle kit; Illumina, San Diego, CA, USA). Using cutadapt software (Martin, 2011) in Python version 2.7, barcode sequences were removed, reads trimmed to a minimum quality score of 20, and then reads less than 250 bp in length were discarded (Code S1). Forward and reverse reads were combined using PANDAseq (Masella et al., 2012) for each sample before all files were merged into one file containing all samples (Code S1). Quantitative Insights Into Microbial Ecology (QIIME) software version 1.9.1 (Caporaso et al., 2010) in Python version 2.7 was used to analyse the data. Operational Taxonomic Units were assigned with 97% clustering to the Greengenes database version 13.8 (DeSantis et al., 2006) for 16S rRNA gene alignment. Sequences aligning to mitochondria or chloroplast sequences were screened for and removed from the dataset. Custom scripts in QIIME were used to analyse relative taxonomic abundance, and alpha and beta diversity measures (Code S1) at a sequencing depth of 10,000 reads per sample. Pairwise comparisons of beta diversity measures (weighted (Lozupone et al., 2007) and unweighted (Lozupone & Knight, 2005) UniFrac) were made using 999 Monte Carlo permutations (MCP). The linear discriminant effect size (LEfSe) (Segata et al., 2011) algorithm was performed to identify taxonomic groups associated with the variables measured (p < 0.01, linear discriminant analysis (LDA) score (log10 ≥ 2)). To be included in the results, each variable must have met the inclusion criteria (p < 0.01, LDA score (log10 ≥ 2)) within each child, as well as when averaged across all three children. Higher taxonomic levels were excluded where it was assumed that a lower taxonomic level was accountable for the observed change. These situations were where a higher taxonomic level had a less significant or equal change in relative abundance compared to a lower taxonomic level classified to the higher taxonomy by LEfSe analysis. However, if the higher taxonomic level had a more significant p-value it was retained. Kruskal–Wallis tests to compare read counts (significantly different if p < 0.05) were performed in R version 3.4.2 (R Core Team, 2017) and graphs were generated using the ggplot2 package (Wickham et al., 2018). All data are provided as Supplemental Information.

Alpha diversity scores (species richness, Shannon and Simpson), generated using standard parameters in QIIME 1.9.1, were analysed by generating linear mixed effect models using the lme4 package (Bates et al., 2018) in R version 3.4.2 (R Core Team, 2017) to identify important predictors of alpha diversity. The lmerTest package (Kuznetsova, Brockho & Christensen, 2017) was used to determine the significance of these model components. Two maximal models were constructed and included all the fixed effects (preservation method, time-to-freezing and stool region) and their interactions with child replicate as a random effect, the first included time as a continuous variable and the second included time as a factor. Backward elimination was used for sequential removal of non-significant variables, to obtain the minimal statistically significant model (Burnham, Anderson & Huyvaert, 2011) (Code S2).

Taxonomic abundance graphs and LEfSe plots generated with QIIME were recreated in R using the ggplot2 package (Wickham et al., 2018). Principal coordinates analysis (PCoA) plots were generated using Emperor Software (Vázquez-Baeza et al., 2013) within QIIME.

Due to the low number (n = 2) of OMNIgene.GUT samples taken, and because OMNIgene.GUT samples were only taken from two out of the three children, these samples were excluded from the above analyses and analysed separately.

Samples processing and microbiome sequencing

In total 87 stool samples were collected for analysis and libraries prepared. After sample exclusion, trimming and alignment (Figs. S1 and S2; Table S1) there was an average of 67,575 (range 19,083–466,807) reads per sample (n = 85).

Microbiome profiles vary between individual children

Each child had a distinct microbiome signature (Figs. 1A and 1B) that was apparent at all taxonomic levels, from phylum (Fig. 2A; Table S2) to genus (Fig. 2B; Data S1) level regardless of the preservation method used and the time-to-freezing duration. The most abundant phyla were Bacteroidetes in child A (40.7%) and child B (36.9%), and Firmicutes in child C (34.1%); followed by Firmicutes in child A (40.1%), Proteobacteria in child B (30.2%) and Bacteroidetes in child C (28.7%). LEfSe identified several bacterial taxa significantly associated with each individual child (Data S2). PCoA analysis using qualitative (presence/absence) differences (unweighted UniFrac, Fig. 1A) confirmed that the clustering of bacterial sequences within individual children was significantly different (MCP for all child comparisons p ≤ 0.001) (Fig. 1A). Children were also significantly different by relative abundance weighted UniFrac (MCP for all child comparisons p ≤ 0.001) (Fig. 1B).

Microbiome profiles vary by stool storage method used

The samples stored as raw stool had a mean average of 80,822 reads per sample (range 26,270–466,807; 35 samples), samples in ethanol had a mean average of 62,983 reads per sample (range 24,215–140,356; 24 samples), and samples in RNAlater had a mean average of 53,981 reads per sample (range 19,083–100,934; 26 samples) (Fig. S3). The number of read counts was not significantly different between preservation methods using a Kruskal–Wallis test.

The within-individual variation between samples stored under different preservation methods was less than that observed between individuals (Figs. 1A, 1B and 2). Intra-storage preservation method microbiota abundance compositions were similar over the time points examined (0–32 h) for each of the preservation methods used (Fig. 2), suggesting relative stability. There were twelve taxonomic groups significantly associated with raw stool (all p < 0.01, LDA score (log10 ≥ 2), Data S3). Eight and eleven taxonomic groups were positively associated with ethanol and RNAlater storage respectively, compared to raw stool alone (all p < 0.01, LDA score (log10 ≥ 2) (Fig. 3; Data S3). Seven of the taxonomic groups that were significantly elevated in relative abundance were shared in the samples stored in ethanol and RNAlater (Fig. 3; Data S3).

Microbiome diversity under different preservation methods were found to differ significantly within each child by qualitative unweighted UniFrac analysis (MCP; all comparisons p ≤ 0.001). The unweighted Unifrac distances within ethanol samples were found to be significantly different compared to within RNAlater samples (child A: p = 0.003, child B: p = 0.041, child C: p = 0.035) or within raw stool samples (child A: p = 0.02, child B: p = 0.038, child C: p ≤ 0.001). UniFrac metrics within raw stool samples were not significantly different to metrics within RNAlater samples apart from in child C (p ≤ 0.001). MCP analysis of unweighted UniFrac comparisons across all three children also revealed significant differences between raw stool and RNAlater storage methods (raw stool vs. raw stool:raw stool vs. RNAlater, p = 0.005; RNAlater vs. RNAlater:RNAlater vs. raw stool, p = 0.01), and raw stool and ethanol storage (raw stool vs. raw stool: raw stool vs. ethanol, p = 0.026) (Fig. 1C) despite distinct separation by child (Fig. 1A).

When UniFrac measures were weighted by relative sequence abundance within each child and as an average of all three children, all preservation method comparisons by MCP were found to be significantly distant from each other (for all comparisons p ≤ 0.001) (Figs. 1B and 1D). Separation within raw stool samples was also found to be significantly different compared to separation within ethanol storage (MCP; child A: p = 0.017, child B: p = 0.041, child C: p ≤ 0.001, all children: p = 0.005). In child C within RNAlater metrics were found to be significantly different by MCP compared to within raw stool metrics (p ≤ 0.001) however, this was not observed in child A or child B. No significant differences were observed by MCP between within RNAlater and within ethanol weighted UniFrac metrics.

Microbiome profiles remain relatively stable over time

Relative bacterial abundance and composition remained relatively stable over time-to-freezing across all storage techniques by LEfSe analysis when time-to-freezing was included as a continuous variable (Fig. 4). However, LEfSe analysis indicated a significantly increased relative abundance of two bacterial genera in raw stool samples at 32 h time-to-freezing when time-to-freezing was included as a categorical variable: Sediminibacterium (p = 0.0016, LDA (log10) = 2.35) (Fig. 5A) and Rummeliibacillus (p = 0.0012, LDA (log10) = 3.25) (Fig. 5B). No significant categorical time-to-freezing effects were identified in ethanol or RNAlater samples by LEfSe analysis. No apparent time clustering was observed by PCoA using UniFrac metrics. No significant differences by weighted or unweighted UniFrac metrics using MCP were observed when comparing 0 h sample metrics to any other time-to-freezing time point or vice versa.

No significant differences in microbiome composition were observed between inner, outer and mixed regions of stool samples

Microbiome profiles of the mixed stool samples were similar to inner and outer stool samples, from phylum through to genus level (Fig. 6). Significant differences between stool regions were not observed in this study by LEfSe regardless of whether the data were analysed by child, storage preservation method or as a whole. No significant associations were generated when LEfSe analysis was performed and no apparent clustering was identified by PCoA analysis of UniFrac diversity metrics. Weighted and unweighted UniFrac comparisons showed no significant differences using MCP.

Modelling indicates storage method influences stool alpha diversity

Linear mixed effect models were constructed to detect variables associated with alpha diversity metrics (Code S2; Table S3). None of the final models identified stool region or time (included as a fixed or a continuous variable) to be significant predictors of alpha diversity. Individual children included as a random effect alone was the only variable shown to influence Shannon diversity and Simpson diversity.

Shannon ~ 1 + (1|child)

Simpson ~ 1 + (1|child)

Preservation method was identified to be a significantly important model component for the model predicting species richness (p = 2.871e−13) with stool stored in RNAlater having the highest average richness, followed by raw stool and stool stored in ethanol (Table S3). Compared to the null species richness model, 50.2% of species richness variation was accounted for by the final species richness model.

Species richness ~ preservation method + (1|child)

OMNIgene.GUT sample performance

Samples stored using OMNIgene.GUT kits had a mean average of 77,089 reads per sample (child B = 102,227; child C = 51,951). The samples clustered by PCoA analysis to the relevant children from which they were taken (Fig. S2). Relative abundance profiles were also representative of the microbiome profiles from each child (Fig. S4).