Different analysis strategies of 16S rRNA gene data from rodent studies generate contrasting views of gut bacterial communities associated with diet, health and obesity

Ethical considerations

This study used 16S rRNA gene sequence data from public datasets published by our research group (Table 1), all of which received approval for the use of the animals.

16S gene sequencing data

We used 16S gene sequencing data from seven of our previous publications dealing with modulation of the gut microbiota using dietary interventions in animal obese models (Table 1). Performing a comparative analysis of data generated by the same research group is advantageous because technical variation is likely to be less compared to the variation obtained from multiple research groups. The sequencing procedure was performed at UC Davis for three studies, while the remaining four studies were sequenced at the Molecular Research LP. The Quick-DNA Fecal/Soil Microbe Miniprep Kit (D6010; Zymo Research) was used in all studies with only slight variations except for the peach study that used the QIAamp^® PowerFecal^® DNA kit (Qiagen). We did not use FastPrep (MP Biomedicals) for bead beating in any of these studies and this is important to clarify because some people consider the use of FastPrep essential for optimal lysis.

Six factors were studied across the samples: ‘dietary treatment’, with 11 levels; ‘study’, with seven levels, one for each study; ‘animal model’, with two levels: mice and rats; ‘sequencing technique’, with two levels: pyrosequencing and MiSeq; ‘obesity status’ at the time of sampling, with two levels: lean and obese; and ‘anatomical site’, with two levels: colon and feces. It is important to keep in mind that in the case that one factor is biologically significantly associated with a different microbiota in nature (e.g., obesity status), the existence of other interactive factors (e.g., dietary treatment) may mask the differences we observe during analysis. Unless otherwise stated, the data was analyzed with QIIME (Caporaso et al., 2010) v.1.8.0. After demultiplexing and quality filtering, we used the sequence file to assign 16S reads to OTUs based on the GreenGenes reference OTU sequence files (v.13.8, August 2013 release) clustered at 97% (99,322 sequences) and 99% similarity (203,452 sequences).

Assignment of 16S sequences to OTUs

Three strategies were used to assign the 16S sequences to OTUs. First, the conservative closed approach that discards sequences for not matching any sequence in the sequence reference data. Second, a de novo approach, which does not use a reference data set to cluster the sequences (Westcott & Schloss, 2015), thus being relatively free to depict the variety of sequences in a sample. Finally, an open approach that combines these two approaches (Rideout et al., 2014), first performing a closed approach followed by clustering of remaining sequences de novo. From this point on, we will refer to these strategies as the name of the clustering method (e.g., closed) and the reference file utilized (e.g., OTUs clustered at 97%). For example, ‘closed97’ and ‘closed99’ will refer to a closed approach using the 97% and the 99% OTU reference files, respectively. We only used the default 97% in the similarity option of the pick_otus.py script but changing this parameter can also drastically affect the results (more on this in “Similarity percentage between 16S rRNA gene sequences” in Supplemental Information).

Chimeras

Chimeras are possible combinations of two or more parent sequences that can inflate the observed OTUs and other diversity parameters, especially for de novo and open-reference strategies to select OTUs. However, there is a lack of consensus in chimera detection and removal and even the existence of true chimeras has been questioned. In this study, we subsampled the pynast alignment of representative sequences (one per OTU) separately for each study and applied the public version of uchime (v.4.2.40) available in Mothur for chimera detection. These analyses were performed on representative sequences from the de novo and open strategies. All the GreenGenes reference OTU sequence files (total: 14 files; OTUs clustered at different similarity percentages, from 99% similarity with 203,451 sequences, to 61% similarity with 22 sequences) were used to explore the behavior of uchime.

Taxonomic and diversity analyses

We combined all OTU tables (one for each of the seven studies) from the closed97 and the closed99 approaches into two separate OTU tables (one for each approach) and used these OTU tables for taxonomic classification and diversity analyses. The relative abundance of each taxon was calculated using all sequences (including possible chimeras and OTUs that appear only once) because it is impossible to detect true chimeras and the possible relevance of low abundant groups in the gut microbiome (Claussen et al., 2017). The unique fraction metric, or UniFrac, is a phylogenetic method for comparing microbial communities based on the phylogenetic distance between sets of taxa in a phylogenetic tree as the fraction of the branch length of the tree that leads to descendants from one environment or the other, but not both (Lozupone & Knight, 2005). Both weighted and unweighted UniFrac distances were used for comparing microbial communities because they can lead to different insights into factors that structure microbial communities (Lozupone et al., 2007). Principal Coordinate Analyses (PCoA) using these UniFrac distances were performed in QIIME and visualized using Emperor. Additionally, we used uniform manifold approximation and projection (UMAP, McInnes et al., 2018), a non-linear dimensionality reduction technique, to confirm the observed clusters, using the umap R package.

Prediction of functional profiling and phenotypes

PICRUSt (phylogenetic investigation of communities by reconstruction of unobserved states, Langille et al., 2013) was used to predict the functional profiling based on the 16S sequences. PICRUSt results were analyzed in STAMP (Parks & Beiko, 2010). Additionally, we used BugBase (Ward et al., 2017) to predict organism-level microbiome phenotypes for each study separately.

QIIME2

QIIME2 was introduced in 2017 on the basis of a plugin architecture that allows third parties to contribute functionality (Bolyen et al., 2019) and has been constantly updated ever since (more than 20 versions of QIIME2 have been released). This paper was conceived and started in QIIME1 in 2018 but because QIIME1 is no longer updated or supported, DADA2 (Callahan et al., 2016) and Deblur (Amir et al., 2017) were also used to select OTUs for a few selected studies in QIIME2, in part based on the observations made by Thompson et al. (2017) in a context of OTUs and sequence variants. One interesting feature of DADA2 is the removal of chimeric sequences.

Statistical analysis

The non-parametric analysis of similarities (ANOSIM) and the Adonis tests were used to determine whether the clustering of samples by a given factor (e.g., study) is statistically significant based on UniFrac distances, using the compare_categories.py script in QIIME with default number of permutations (999). In our experience, these tests usually have low sensitivity (they usually yield low p values even for weak clustering of samples), therefore it is informative to look at both the p values and the percentage of variation explained by the factor. We used the non-parametric Kruskal-Wallis to compare the number of OTUs between the different levels of any given factor, and also to compare the results from BugBase.

Community membership structure

Table 2 shows an accurate numerical impression of relative abundances of all bacterial phyla across different OTU picking strategies, with up to 20% difference in relative abundance of some taxa, depending on the strategy analysis. However, Table 2 lacks a general view of how each taxon is represented across the different levels of each factor. A visual analysis of membership data revealed interesting patterns, for instance, the relative abundances of Firmicutes and Bacteroidetes were equally represented in each level of all six factors studied, while Cyanobacteria was poorly represented in samples from mice (Fig. 1). However, note that these results may also be misleading when considering the very few OTUs that were shared across all samples.

Peach study

The peach study was the only one using 454 pyrosequencing, a different DNA extraction method, and a different primer set. As expected, an increase in quality threshold reduced the number of sequences that passed quality filtering. This is important because it has been suggested that a more stringent quality filtering helps to “reduce the number of spurious OTUs” (Buza et al., 2019), although in practice it is difficult to determine the exact quality threshold to differentiate “true” vs “false” OTUs (Edgar, 2017). We will not discuss this issue for Illumina platforms because defaults have already been established (Bokulich et al., 2013).

Overall the different taxa remained similar in relative abundance but a higher base quality score (qual) threshold of 34 (default: qual 25) had a strong effect on the number of sequences available for OTU picking (62,971 sequences with qual 25 vs. 11,553 with qual 34) and consequently on the number of the OTUs discovered (758 OTUs with qual 25 vs 442 OTUs with qual 34 in the closed97 approach). The difference in qual threshold also had an effect in the proportions of Bacteroidetes that went down from ∼37% (qual 25) to 15% (qual 34), and Firmicutes that went up from ∼50% to 70%. This discrepancy was not related to the lower number of sequences available for OTU picking because a lower rarefaction in all other analyses (e.g., with qual 25) revealed similar relative proportions compared to the analysis with higher rarefaction depth. A higher base quality score threshold also affected the presence of some low abundant groups (e.g., at qual 34 Deferribacteres and Fusobacteria were not detected using all approaches).

In the peach study, Bacteroidetes displayed the highest standard deviation (SD) and showed the biggest difference (∼5%) in relative abundances, particularly between the open97 and the open99 approaches (Table 2). The difference between the lowest and the highest value was minimal for Firmicutes (∼2%), Proteobacteria (∼2%) and others (Table 2). Tenericutes showed the highest SD/average ratio (68.7), which implies that the variability was proportionally higher in this taxon, and Firmicutes the lowest (1.7), which implies that the variability was proportionally lower in this taxon. The detected phyla varied from 10 (closed97) to 13 (de novo97 and open97).

Wheat study

In the wheat study, Firmicutes displayed the highest SD and showed the biggest difference (∼20%) in relative abundances, particularly between the closed97 and the closed99 approaches (Table 2). This contrasts heavily with the biggest difference of ∼5% observed in the peach study. This difference in closed approaches was also noticeable in Bacteroidetes (∼11% difference), Proteobacteria (∼10% difference), and others (Table 2) and this was not related to rarefaction depth. This difference was not noticeable for the open and the de novo approaches where Firmicutes, Bacteroidetes, Proteobacteria and others showed very similar proportions in all cases (Table 2). Verrucomicrobia showed the highest SD/average ratio (107.9) and Firmicutes the lowest (11.3). The number of detected phyla varied between 10 (closed97) and 14 (de novo).

Quinoa study

Similar to the wheat study, in the quinoa study Firmicutes displayed the highest SD and showed the biggest difference (∼15%) in relative abundances, particularly in the closed approaches (Table 2). This difference in closed approaches was also noticeable in Bacteroidetes (∼11% difference) and minimal for Proteobacteria (∼2% difference) and others, and this was again not related to rarefaction depth. Also similar to the wheat study, this difference was not noticeable for the open and the de novo approaches where Firmicutes, Bacteroidetes, Proteobacteria and others showed very similar proportions in all cases (Table 2). Verrucomicrobia showed the highest SD/average ratio (101.8) and Firmicutes the lowest (7.3). The number of detected phyla varied from 8 (closed97) to 12 (de novo).

Barley study

In the barley study, Bacteroidetes displayed the highest SD but Firmicutes showed the biggest difference (∼9%) in relative abundances, particularly in the closed approach (Table 2). This difference in closed approaches was minimal in Bacteroidetes (∼3% difference), Proteobacteria (∼2% difference) and others and this was again not related to rarefaction depth. Also similar to the other studies, this difference was not noticeable in the de novo and open approaches (Table 2). Verrucomicrobia showed the highest SD/average ratio (91.9) and Firmicutes the lowest (5.5). The number of detected phyla varied from 8 (closed97) to 12 (de novo).

Cherry study

In contrast to the wheat, quinoa and barley studies, where the closed97 and the closed99 approaches showed different abundance of taxa, in the cherry study there was good agreement between the closed97 and the closed99 approaches for Firmicutes, Bacteroidetes, Proteobacteria and other groups, but these approaches showed dissimilar proportions compared to the open and the de novo approaches (Table 2). The cherry study was also interesting because the de novo approach showed the presence of bacterial groups (e.g., Nitrospirae, Chlorobi, Planctomycetes) that were not detectable using the closed and the open approaches, and this was likely not related to rarefaction depth because the open approach used similar thresholds. Bacteroidetes displayed the highest SD and showed the biggest difference (∼7%) in relative abundances. Tenericutes showed the highest SD/average ratio (100.9) and Firmicutes the lowest (3.4). The number of detected phyla varied from 14 (open99) to 22 (de novo). The use of DADA2 and Deblur in QIIME2 showed 1,329 and 1,263 OTUs, respectively. This is about half the lower number of OTUs (2,439 with closed97) obtained with the other approaches (Table 2).

Raspberry study

Unlike the other studies discussed above, the raspberry study was very interesting because it had the lowest variation among the different OTU picking approaches (0.56% difference in Bacteroidetes and 0.55% difference in Firmicutes between the lowest and the highest results), which was also reflected in the SD/average ratio (highest: 32.8 for Cyanobacteria, lowest: 1.2 for Firmicutes) (Table 2). Moreover, both the open97 and the open99 approaches detected as many as 41 different taxa at the phylum level, while the number of taxa in all other studies using the same approach only varied from 8 to 22. This strongly suggests that the microbiota data contained within each study is unique and may sometimes contain a high proportion of taxa that go unnoticed in other studies. This is very important in a context of the role of rare taxa in maintaining the stability of ecosystems (Jousset et al., 2017). DADA2 and Deblur showed 791 and 721 OTUs, respectively. This contrast with the numbers (2,751-92,486 OTUs) obtained from all approaches (Table 2).

Apple study

In the apple study, Bacteroidetes showed the highest variation and also the biggest difference (∼11%) especially between the closed and the de novo approaches (Table 2). This difference between the two approaches was also noticeable in Firmicutes and Tenericutes (Table 2). Tenericutes showed the highest SD/average ratio (77.3) and Proteobacteria the lowest (4.9). The number of detected phyla varied from 15 to 35, thus making the apple study the second study more variable after the raspberry study. The use of DADA2 and Deblur in QIIME2 showed 1,357 and 1,149 OTUs, respectively. This contrast with the numbers (2,095-153,681 OTUs) obtained from all approaches (Table 2).

Firmicutes/Bacteroidetes ratio

The first paper that discussed about the possible usefulness of the Firmicutes/Bacteroidetes ratio was published by Ley et al. (2005), where the showed that obese mice had 50% more Firmicutes, with a proportionally lower abundance of Bacteroidetes, thus leading to a higher Firmicutes/Bacteroidetes ratio compared to lean mice. However, each phylum is composed by hundreds of different species, and therefore the ratio between the two has little significance as discussed elsewhere (Delzenne & Cani, 2011). In this current study, this ratio varied from 1.1 to 3.3 due to the differences in relative abundance in the two phyla (Table 2).

The phylum Verrucomicrobia

Firmicutes and Bacteroidetes are usually the most abundant phyla in the gut microbiota and they have attracted most of the attention. However, other taxa deserve attention to better comprehend the functioning of the gut microbial ecosystem. A. muciniphila is a taxon that has generated interest as a new generation probiotic candidate to help obese patients and is considered to be a member of Verrucomicrobia based on 16S gene analysis. Akkermansia was detected in ∼15% the samples (25/164) with an average of 3% and it was more represented in colon samples (Fig. 1). Interestingly, the whole Verrucomicrobia phylum also had an average of 3%, which implies that most or all Verrucomicrobia was represented by Akkermansia. If we translate the proportions of 16S reads into numbers, a 3% would represent approximately 300 million cells (3 ×10⁸) in a hypothetical environment of 1 ×10¹⁰ cells/g of intestinal contents. While this conversion is not necessarily accurate, the high numbers may bear some ecological relevance if one considers that as low as 1,000 cells of other microbes are enough to thrive (Nilsson, Kari & Steele-Mortimer, 2019). In three studies, the abundance of Verrucomicrobia was >5 times higher using a closed approach compared to the other approaches (Table 2).

The phylum Cyanobacteria

Another group of interest for gut health is Cyanobacteria (Barcena et al., 2019), which is often considered a contaminant and removed from 16S gene analysis. However, Ley et al. (2005) showed a deep-branching clade of Cyanobacteria in the guts of animals and mentioned that they may represent descendants of non-photosynthetic ancestral Cyanobacteria that have become adapted to life in the mammalian gut. Here, Cyanobacteria was detected in ∼50% the samples (74/164) and was always low in relative abundance (<0.5%) (Table 2). While 0.5% may be considered low, this percentage implies approximately 50 million (5 ×10⁷) cells in the same hypothetical environment of 1 ×10¹⁰ cells/g of intestinal contents as discussed above.

Analyses of UniFrac distances

UniFrac analyses from the closed97 approach

The first two coordinates explained 18% of the variation using unweighted distances (Fig. 2), which is similar to other studies that have shown that the first two coordinates explained 23% of the variation in a data set of unweighted UniFrac distances from samples of the human microbiota (Lozupone et al., 2013). Accordingly to ANOSIM tests on unweighted distances, study, sequencing procedure, animal model, and dietary treatment were the factors most highly associated with the differences in the bacterial communities (Fig. 2, Table 4). UMAP confirmed the clustering by animal model and study (Supplemental Information). This is interesting because others have also shown a strong study effect (Lozupone et al., 2013).

Table 4:

Summary of results for all samples (n = 162) from the Adonis and ANOSIM tests for comparing categories using UniFrac data from the closed97 approach.

	Adonis		ANOSIM
	Unweighted	Weighted	Unweighted	Weighted
Study	P < 0.001	P < 0.001	P = 0.001	P = 0.001
	R2 = 24.8%	R2 = 19.8%	R = 66.3	R = 31.5
Model	P < 0.001	P < 0.001	P = 0.001	P = 0.145
	R2 = 8.8%	R2 = 7.4%	R = 56.7	R = 4.5
Sequencing	P < 0.001	P < 0.01	P = 0.001	P = 0.612
	R2 = 5.5%	R2 = 3.2%	R = 57.8	R =  − 3.2
Treatment	P < 0.001	P < 0.001	P = 0.001	P = 0.001
	R2 = 20.4%	R2 = 21.1%	R = 34.2	R = 20.1
Obesity	P < 0.001	P < 0.001	P = 0.023	P = 0.001
	R2 = 1.9%	R2 = 4.4%	R = 9.5	R = 23.6
Site	P < 0.001	P < 0.01	P = 0.702	P = 0.868
	R2 = 3.8%	R2 = 2.2%	R =  − 1.9	R =  − 4.7

DOI: 10.7717/peerj.10372/table-4

Notes:

According to QIIME documentation, the R² value (effect size) calculated with Adonis test shows the percentage of variation explained by the factor (e.g., the factor ‘study’ explained 24.8% of the variation in the unweighted distances), and the R statistic calculated with ANOSIM test reflects the degree of dissimilarity (an R value near 1, or 100, means that there is dissimilarity between the groups, while an R value near 0 indicates no significant dissimilarity between the groups). A rarefaction depth of 100 sequences per sample to account for as many samples as possible (only two samples were left out using this rarefaction depth). A total of 999 permutations were used to calculate the statistics.

In this study, the first 2 coordinates explained 34% of the variation using the weighted UniFrac distances, which is about twice the variation explained by the first two axes using unweighted UniFrac, but produced a very similar clustering of samples. Accordingly to ANOSIM tests, the effect of study, animal model, and sequencing procedure explained much smaller proportion of the variability in the data using weighted distances, while obesity status explained a higher proportion of the variation (23.6% vs 9.5% using unweighted distances, Table 4).

UniFrac analyses from the closed99 approach

The analysis of unweighted UniFrac distances from the closed99 approach revealed similar results compared to the closed97 regarding the relative contribution of each factor in explaining the variability in the data (Fig. 3 and Table 5). The first two coordinates explained 25% of the variation (Fig. 3), which is even closer to the results of Lozupone et al. (2013). Interestingly, Akkermansia was not part of the ten most abundant groups (Fig. 3) and the factor ‘study’ explained most of the dissimilarity in the communities using the unweighted UniFrac data accordingly to the ANOSIM test (Table 5). In our experience it is uncommon that one single factor explains that much of the variability between microbial communities. Similar to the analysis using the closed97 approach, the clustering of samples was very similar but the effect of ‘study’, ‘animal model’ and ‘sequencing procedure’ at explaining the variability in the data was lower using weighted UniFrac, while obesity status explained much more of the variability in the data (Table 5).

Table 5:

Summary of results for all samples (n = 163) from the Adonis and ANOSIM tests for comparing categories using UniFrac data from the closed99 approach.

	Adonis		ANOSIM
	Unweighted	Weighted	Unweighted	Weighted
Study	P < 0.001	P < 0.001	P = 0.001	P = 0.001
	R2 = 33.2%	R2 = 33.4%	R = 80.1	R = 43.5
Model	P < 0.01	P < 0.001	P = 0.001	P = 0.016
	R2 = 11.9%	R2 = 10.4%	R = 53.9	R = 9.8
Sequencing	P < 0.001	P < 0.001	P = 0.001	P = 0.529
	R2 = 6.7%	R2 = 3.7%	R = 61.6	R =  − 2
Treatment	P < 0.001	P < 0.001	P = 0.001	P = 0.001
	R2 = 24.6%	R2 = 28.2%	R = 33.8	R = 25
Obesity	P < 0.001	P < 0.001	P = 0.237	P = 0.001
	R2 = 1.8%	R2 = 5.5%	R = 3.1	R = 20.6
Site	P < 0.001	P < 0.01	P = 0.960	P = 0.953
	R2 = 5.1%	R2 = 2.4%	R =  − 5.3	R =  − 6.3

DOI: 10.7717/peerj.10372/table-5

Notes:

According to QIIME documentation, the R² value (effect size) calculated with Adonis test shows the percentage of variation explained by the factor (e.g., the factor ‘study’ explained 33.2% of the variation in the unweighted distances), and the R statistic calculated with ANOSIM test reflects the degree of dissimilarity (an R value near 1, or 100, means that there is dissimilarity between the groups, while an R value near 0 indicates no significant dissimilarity between the groups). A rarefaction depth of 440 sequences per sample to account for as many samples as possible (only one sample was left out using this rarefaction depth). A total of 999 permutations were used to calculate the statistics.

Prediction of functional profiling

PICRUSt revealed a total of 217 features that were significantly different between different studies, 211 between different treatments, 145 between sequencing techniques, 136 between mice and rats, 69 between lean and obese, and 32 between feces and colon contents (adjusted P < 0.05, see ‘PICRUSt results’ in Supplemental Information). This is interesting because the factor ‘study’ was also associated with the highest amount of variability. The weighted Nearest Sequenced Taxon Index (weighted NSTI) score was developed to summarize the extent to which microorganisms in a given sample are related to sequenced genomes (e.g., a NSTI score of 0.03 indicates that the genomes of the microbes were well represented and that the average microbe in the sample can be predicted using a relative from the same “species”). Higher NSTI values indicate a lower accuracy of the predictions. Here, NSTI scores varied widely, from a minimum of 0.03 (sample W11, feces from lean mice sequenced with MiSeq, wheat study) to a maximum of 0.34 (sample W32, feces from obese mice on quinoa sequences with MiSeq, quinoa study). There was no indication to suggest that any group of sequences were associated with lower or higher NSTI scores.

Prediction of organism-level microbiome phenotypes

BugBase revealed interesting differences in phenotypes within each study. For example, obese rats in a high-fat diet had higher proportions of bacteria with the potential of forming biofilms among the different treatments in the apple study, lean mice had higher proportions of bacteria containing mobile elements in the barley study, and obese mice had higher proportions of stress-tolerant bacteria in the cherry study (more on this in ‘BugBase results’ in Supplemental Information). Overall, these results also indicate that each study showed unique peculiarities, ultimately derived from the molecular composition of the 16S gene sequences at any given time-point within each particular study.