Influence of fecal collection conditions and 16S rRNA gene sequencing protocols at two centers on human gut microbiota analysis

Background To optimise fecal sampling and analysis yielding reproducible microbiome data, and gain further insight into sources of its variation, we compared different collection conditions and 16S rRNA gene sequencing protocols in two centers. Fecal samples were collected on three sequential days from six healthy adults and placed in commercial collection tubes (OMNIgeneGut OMR-200) at room temperature or in sterile 5 ml screw-top tubes in a home fridge or home freezer for 6-24 h, before transfer at 4°C to the laboratory and storage at - 80°C within 24 hours. Replicate samples were shipped on dry ice to centers in Australia and the USA for DNA extraction and sequencing of the V4 region of the 16S rRNA gene, using different PCR protocols. Sequences were analysed with the QIIME pipeline and Greengenes database at the Australian center and with an in-house pipeline and SILVA database at the USA center. Results Variation in gut microbiome composition and diversity was dominated by differences between individuals. Minor differences in the abundance of taxa were found between collection-processing methods and day of collection. Larger differences were evident between the two centers, including in the relative abundances of genus Akkermansia, in phylum Verrucomicrobiales, and Bifidobacteria in Actinobacteria. Conclusions Collection with storage and transport at 4°C within 24 h is adequate for 16S rRNA analysis of the gut microbiome. However, variation between sequencing centers suggests that cohort samples should be sequenced by the same method in one center. Differences in handling, shipping and methods of PCR gene amplification and sequence analysis in different centers introduce variation in ways that are not fully understood. These findings are particularly relevant as microbiome studies shift towards larger population-based and multicenter studies.

particularly relevant as microbiome studies shift towards larger population-based and multicenter studies.

Background
The fecal or 'gut' microbiome is shaped strongly by diet but also by the host genotype, age, hygiene and antibiotic exposure, and is altered in many pathophysiological states [1], [2]. The composition of gut microbiota differs greatly between individuals [3] and therefore maximizing the detection of biological and disease-related changes requires minimization of variation due to methods of sample collection and analysis.
Previous studies have shown that fecal microbial composition overall was not altered when DNA was extracted from a fresh fecal sample compared to a sample that had been immediately frozen and stored at -80 °C for up to 6 months [4,5]. Storage at different temperatures for varying times has been compared with immediate freezing and storage at -80 °C. One study [5] reported a decrease in Bacteroidetes and an increase in Firmicutes phyla after 30 minutes at room temperature, but the majority [4,[6][7][8][9][10][11] have found that storage at room temperature for at least a day, or at 4, -20 or -80 °C for up to 14 days, had little effect on the relative abundance of taxa. Moreover, microbial composition was not significantly altered in DNA extracted from fecal occult blood test cards that had been at room temperature for three days [12]. Recent studies have also evaluated the OMNIgene®•GUT (OMR-200) collection and liquid storage tube, which is reported to stabilize DNA at room temperature for 14 days [13]. Samples immediately collected into these tubes and stored for three days at room temperature showed little difference in microbial composition by 16S rRNA gene sequencing compared with samples immediately frozen at minus 80 °C [7]. A similar result was obtained when samples in the tubes, stored for 1-28 days at room temperature, were compared with corresponding samples from which DNA had been freshly extracted [14,15]. However, the relative abundance of Bacteroides increased after seven days and in infants (with lower microbial diversity than adults) significant divergence from fresh samples was observed after 14 days storage.
As microbiome studies expand into larger populations at multiple sites, stringent quality control remains critical. With this in mind and to optimise analysis of the gut microbiome in a multi-site, longitudinal pregnancy-birth cohort study [16], we evaluated different collectionprocessing methods on three sequential daily fecal samples from six individuals, and 16S rRNA gene sequencing in two centers.

Methods
Six healthy adult volunteers, three males and three females aged 35-70, provided fecal samples on three successive days. Multiple aliquots were taken from each bulk sample and those for bacterial microbiome analysis were stored by one of four methods A-D ( Figure 1).

Collection-processing methods
Method A: individuals placed aliquots of feces into 6 x OMNIgene®•GUT (OMR-200) [13] tubes as per manufacturer's protocol. These were stored at room temperature for 6-24 h before delivery to the laboratory for transfer to sterile 5 mL screw cap tubes and storage at minus 80°C. Method B: individuals placed aliquots of feces into 6 x sterile 5 ml screw cap tubes, which were stored in the home freezer for 6-24 h before delivery to the laboratory in an insulated container for storage at minus 80°C. Method C: individuals placed aliquots of feces into 6 x sterile 5 m screw cap tubes, which were stored in the home refrigerator for 6-24 h before delivery to the laboratory in an insulated container for storage at minus 80°C. Method D: individuals placed a bulk fecal sample into a sterile 70 ml collection jar, which was stored in home refrigerator for 6-24 h before delivery to the laboratory in an insulated container, transfer of aliquots into 6 x sterile 5 m screw cap tubes and storage at -80°C.
In the laboratory, samples were handled under sterile conditions in a Biosafety Level 2 cabinet.
This collection-processing procedure yielded a total of 432 samples (24 from each of 6 individuals per day for 3 days) ( Figure 1). After 2-4 weeks at -80°C, three sample aliquots per person-day-method (n=216) were transported on dry ice to sequencing centers in Australia

DNA sequencing methods
Samples were thawed on ice and DNA extracted at both WEHI and BCM with the MoBio PowerSoil kit (MoBio Laboratories, Carlsbad, CA), as used in the Human Microbiome Project [17]. At WEHI, the V4 hypervariable region of the bacterial 16S rRNA marker gene (16Sv4) was PCR-amplified in duplicate with primers 515F-OH1 and 806R-OH2. Analogues of these are described, respectively, by the common name U515F, new name S-*-Univ-0515-a-S-19, and the common name 806R, new name S-D-Bact-0787-b-A-20 in [18]. These primers (GTGACCTATGAACTCAGGAGTCGGACTACNVGGGTWTCTAAT) and  [20,21]. The primers used for amplification contained adapters for MiSeq sequencing and single-end barcodes allowing pooling and direct sequencing of PCR products [21]. PCR-amplified amplicons were normalized by concentration before pooling and sequences were generated in one lane of a MiSeq instrument using the v2 kit (2x250 bp paired-end protocol). The 216 samples were analysed in a single pool.

Bioinformatics methods
Two different bioinformatics pipelines were applied to the sequence data, Pipeline W at WEHI and Pipeline B at BCM. In Pipeline W, sequences were clustered into operational taxonomic units (OTUs) with 97% similarity using QIIME (Version 1.8.1) [22][23][24][25], and taxonomically classified by aligning the representative sequences to the Greengenes 13_08 database [26].
Paired-end sequences were assembled with PEAR [27] with parameters -v 100 -m 600 -n 80, where -v is minimum overlap, -m is maximum assembled length and -n is minimum assembled length. WEHI index sequences were extracted with QIIME script extract_barcodes.py, and bases up to and including the 533F to 805R V4 region amplicon primers were trimmed. BCM sequences were supplied as trimmed sequences with a separate barcode index file. QIIME's split_libraries_fastq.py script was used for quality filtering, with phred quality scores required to be above 29, and 90% of a read's length required to have consecutive, high-quality base calls. In order to minimise differences between WEHI and BCM sequences, the WEHI sequences were further filtered by aligning to the SILVA 123 16S rRNA gene database using MOTHUR v1.38.1 [28,29] with start=8390 and end=17068, and removing sequences that were not in this position.
QIIME's open-reference OTU-picking was used with the Uclust algorithm [22,23] to form clusters at 97% similarity. The representative sequences from each OTU were aligned with gaps to a reference set using QIIME's implementation of PyNAST, then filtered for chimeric sequences using UChime with default settings [30].
After making filtered OTU In Pipeline B, read pairs were de-multiplexed based on the unique molecular barcodes using Cassava 1.8.4 from Illumina, and reads were merged using USEARCH v7.0.1090 [23], allowing zero mismatches and a minimum overlap of 50 bases. Merged reads were also trimmed at first base with Q5. In addition, USEARCH was also used to apply a quality filter to the resulting merged reads and reads containing above 0.05 expected errors were discarded.
16S rRNA gene sequences were clustered into OTUs at a similarity cutoff value of 97% using the UPARSE algorithm [33]. OTUs were mapped to an optimized version of the SILVA Database [18] containing only the 16Sv4 region to determine taxonomies. A rarefied OTU table was constructed from the output files generated in the previous two steps for downstream analyses of alpha diversity, beta diversity and phylogenetic trends [34]. Pipeline B was applied only to BCM data.

Results
For WW protocol, library sizes after quality filtering, clustering, and combining PCR replicates ranged from 30,000 to 250,000 sequences per sample, with a median of 67,000 ( Figure S1A); sequences clustered into 12,652 OTUs of minimum size 20. For WB, library sizes ranged from 5,000 to 56,000, with a median of 27,000 ( Figure S1B); sequences clustered into 3,675 OTUs of minimum size 20. The differences between library sizes reflect the differences between technical replicates (three for WEHI, one for BCM), as well as sequencing and filtering differences.

Taxonomic overview
At the phylum level, samples were dominated, as expected, by Bacteroidetes and Firmicutes.
The mean summed proportion of these two phyla was 94%, varying from 71.9% to 99.7% between individual samples. Likewise, a single order from each of these two phyla,

Effect of collection-processing method on taxonomic analysis
With the WW protocol, testing for differential abundances between collection-processing methods revealed no significant differences in counts by phylum, order or family (Table 1, Figure 3A). The differences were tested using DESeq2 with a design controlling for the effect of person and day. A very small number of OTUs (0.04%) were different under collectionprocessing Method A. Table 1. Differences by collection-processing method at three taxonomic levels (analysis WW). accounted for only 2%. Overall, alpha diversity was slightly lower from Method A ( Table 2).   Fig 2B). (C) Mean log (standardised count) plotted against the mean over the collection-processing methods, and a linear regression applied. Method A has the greatest average deviation from the linear model.

Effect of collection-processing method on library size
Collection-processing methods were compared after quality filtering, barcode extraction and clustering. Methods did not differ significantly in the number of sequences identified by the WW protocol. Batch effects were more significant (p<10 -5 ) than collection-processing method, but batch and method together contributed less than 5% of the variation in library size. In summary, different collection-processing methods for fecal samples were associated with only minor differences in the microbiota predicted by 16S rRNA gene sequencing. Method A with the OMNIgeneGut OMR-200 device had minor taxonomic differences compared with Methods B-D that involved immediate freezing or refrigeration. The number of DNA sequences extracted per sample was not different by collection-processing method, and no major differences in abundance were apparent at the phylum level. In the WW analysis, no method was associated with significantly more distance between samples at the OTU level. Thus, different collection-processing methods were associated with minor variations but no consistent bias. Overall, the differences between individuals were much larger than those introduced by the different collection and processing methods.

Effect of sequencing center
Sequences generated at WEHI and BCM were analyzed with pipeline W. The most abundant phyla were similar, but the proportions of less abundant phyla and higher taxonomic resolution differed. For example, the mean proportion of genus Akkermansia in the order Verrucomicrobiales was greater in WB (0.7%) than WW (0.02%). The proportion of Bacteroides was lower in some individual samples for BCM than WEHI ( Figure 4A, also S2).    Initial bioinformatics analysis was performed separately on the WEHI and BCM datasets. For better comparison of the taxonomies, pipeline W was re-applied to a data set comprising the BCM sequences and one of the three WEHI technical replicates ( Figure 5). The ordination plot shows 'batch' effects between the two sequencing centers, and greater between-sample differences in the BCM data. DESeq2 was used to make generalized linear models for the counts at phylum, order and OTU levels ( Table 3). The model included individual ID, day and collection-processing method as factors. At the phylum level, the largest change was in the Verrucomicrobia. At the OTU level, 3% of OTUs were significantly different ( Figure S4, Additional data S1). Most of the differentially abundant OTUs belonged to the orders Clostridiales (63%) and Bacteroidales (31%). The direction of change in OTUs was not consistent, and there were no significant differences in counts for Clostridiales and Bacteroidales between WEHI and BCM. With WB, collection-processing Methods A and B were taxonomically different, with a decrease in Actinobacteria in method A and an increase in Lentisphaerae (although counts were very low) (p < 0.001, Additional Table S1). Lentisphaerae were also increased in method A compared with methods C and D (p < 0.05). There were significant differences between Bray-Curtis distances between samples in WB data (p < 0.001), with collection-processing method A associated with smaller differences between samples from the same individual than methods B, C and D (Additional Table S2). Collection-processing method D resulted in fewer sequences than A or C (p=0.01) with the WB protocol, but the difference was small compared with total variation. ( Figure S5).

Effect of bioinformatics pipeline
BCM sequence data were also analysed with Pipeline B. This used the same fastq files of sequenced DNA of the 16Sv4 region as in WB, but allocated OTUs and assigned taxonomies according to a BCM protocol incorporating a version of the SILVA database [18]. Library sizes varied from 2,300 to 23,000 sequences per sample with a median of 14,000. Sequences were clustered into 463 OTUs, fewer than generated by Pipeline W. Pipeline W used the UCLUST algorithm for clustering, which is known to produce more OTUs [35]. The differences in OTU formation and classification resulted in similar descriptions to protocol WB at the phylum level, but with some differences at the order level ( Figure S6A). The UniFrac distances clustered the samples in a very similar manner to protocol WB ( Figure S6B). The effect of the collection-processing method was stronger than for WW ( Figure S6C), and was similar to WB (Additional Table S3). Differential abundance analysis with DESeq2 indicated a difference in Actinobacteria counts between collection-processing methods, with method A giving lower proportions of Actinobacteria than B or C (p < 0.01). The different bioinformatic approaches gave very different OTU numbers, and some differences in microbiota designations due to differences between the reference databases, for example regarding the status of Akkermansia, and the polyphyletic Clostridia. Under protocol WB, some effects of collection-processing methods were apparent.

Collection
Actinobacteria proportions were reduced in Method A samples compared with the other collection-storage methods, and beta-diversity was also reduced. The reduction in distances between samples for an individual with collection-processing method A possibly indicates that the OMR-200 device mitigated the effects of storage and transport.
The differences between sequencing centers could be due to one or more of several factors.
Although the samples were transported and arrived on dry ice, covert effects of shipping conditions can't be excluded. Similarly, differences in sample handling at the laboratory level can't be excluded. One clear difference between the centers was in the sequencing methods.
Although both centers used the same primers for the primary PCR of the V4 amplicon of the 16S rRNA gene, the approaches then diverged. BCM barcoded the reverse primer for each sample and ran a one-step PCR with 32 cycles of amplification, whereas WEHI added overhang sequences to both primers to subsequently introduce both Illumina sequencing adaptors and dual index barcodes for each amplicon, and ran a two-step PCR with 45-cycles of amplification.
In addition, BCM sequenced a single library pool whereas WEHI sequenced triplicated libraries in separate batches on three different days. These differences in methodology are likely to have contributed to inter-center differences in taxonomic composition.

Conclusions
Collection-processing methods and day of collection contributed to only minor variation in fecal microbiome composition and diversity, the major variation being at the level of the individual. However variation, including at the phylum level, was evident between the two sequencing centers and is likely to be related at least in part to difference in PCR design and conditions. Collection with storage and transport at 4°C within 24 h is adequate for analysis of the gut microbiome, but variation between sequencing centers indicates that cohort samples should be sequenced by the same methods in one center. These findings are relevant to the quality control of microbiome studies, in particular to larger, population-based multi-site studies.

Figure S1
Number

Figure S4
Scatter plot of OTUs showing log 2 fold changes of W counts over B counts versus the mean normalized count. Red dots are OTUs with adjusted p < 0.05, triangles are points with fold changes outside the y-axis limits. The blue lines illustrate the null hypothesis that any change is less than 2-fold. Although the numbers of significant changes are evenly split between increase and decrease, W counts tend to be increased in small OTUs and B counts in larger OTUs.    Bars are colour-coded by phyla using the same colours as in Figure 2. Use of the Silva database for taxonomic assignment has introduced genus labels Prevotellaceae_UCG_001 and Lachnospiraceae_UCG_008. (B) Beta diversity using UniFrac distances between samples. Axes have been reflected to give approximately the same orientation of clusters as for Figure 2C. Directions in NMDS are arbitrary, so the positioning and rotation of clusters does not indicate a real change in the UniFrac distances.
(C) Log of standardised counts (scaled by library size) of four phyla for the four methods for each individual, equivalent to Figure 3A. Points show mean and bars show standard deviation for each individual and collection-processing method (n=9).