3.1 Sampling Location Yields a Variability in Sequencing Depth
The raw data from sequencing generated a total of 560,935 reads, with an average of 20,033 reads per sample. Total read depth per sample was limited to an arbitrary minimum of 7,435 to ensure adequate read depth in any given sample [50], and thus was the cutoff for inclusion in further analysis. Out of the 40 libraries sequenced, we excluded six libraries as the to the read threshold. To retain the paired nature of our analysis, we, therefore, proceeded with fourteen paired samples for which both cloacal swab and cecal data was found. A summary of the complete information can be found in the supplementary tables 1 and 2. The cecal samples had an average of 22,968 reads (IQR 20,524 − 24,516 reads), whereas the cloacal swabs had an average of 17,099 reads (IQR 9,946 − 22,378 reads). The average Good’s coverage was 99.64%, (SD 0.315%), showing that the retained datasets had adequate sequence coverage to sample OTUs. Summary statistics for both the cecal and cloacal datasets are given in Table 1.
Table 1
Descriptive Statistic | Cecal Content Reads | Cecal Content Good’s Coverage | Cloacal Swab Reads | Cloacal Swab Good’s Coverage |
Sample Size | 14 | 14 | 14 | 14 |
Minimum Reads | 18428 | 99.63% | 7435 | 98.98% |
1st Quartile | 20524 | 99.71% | 9946 | 99.34% |
Median | 22926 | 99.74% | 16608 | 99.66% |
Mean | 22968 | 99.73% | 17099 | 99.55% |
3rd Quartile | 24516 | 99.74% | 22378 | 99.77% |
Maximum | 29053 | 99.85% | 31859 | 99.86% |
IQR | 3992 | 0.037% | 12432 | 0.434% |
Range | 10525 | 0.21% | 24424 | 0.88% |
Table 1. Descriptive statistics of the sequencing data for the fourteen cecal content and fourteen cloacal swabs. Column titles are in bold font.
3.2 Broad differences between cecal and cloacal microbiota members
Overall, the twenty-eight samples yielded 1603 OTU’s assigned to 82 Families. The top three families based on relative abundance were Peptostreptococcaceae (12.623%, Phylum Firmicutes), Ruminococcaceae (10.802%, Phylum Firmicutes), and Lactobacillaceae (8.909%, Phylum Firmicutes). Twenty-two families were represented in the cloacal swab samples with eighteen families represented in the cecal content samples, as shown in Figs. 1 and 2, respectively. There were eleven families shared cecal and cloacal samples, with seven families (Atopobiaceae, Bacteria_unclassified, Bifidobacteriaceae, Clostridiales_unclassified, Clostridiales_vadinBB60_group (OTU 0045), Gastranaerophilales_fa (OTU 0010), and Helicobacteraceae (OTU 0004)) unique to cecal content and eleven families (Actinomycetaceae, Clostridiales_vadinBB60_group (OTU 0047), Corynebacteriaceae, Enterobacteriaceae, Enterococcaceae, Gastranaerophilales_fa (OTU 0013), Mollicutes_RF39_fa, Pasteurellaceae, Peptostreptococcaceae, and Planococcaceae) unique to cloacal swabs.
The fourteen cecal content samples yielded 1424 total OTU’s from 54 total Families with Ruminococcaceae (23.062%, Phylum Firmicutes), Barnesiellaceae (12.700%, Phylum Bacteroidetes), and Rikenellaceae (7.841%, Phylum Bacteroidetes) as the top three most abundant families in the cecal content samples. The fourteen cloacal swab samples yielded 905 total OTU's from 79 total Families. The top three most abundant families in the cloacal swab samples do not seem uniform with regards to the distribution of families, nor is there an observable pattern. The top three most abundant families in the cloacal samples were Peptostreptococcaceae (12.912%, Phylum Firmicutes), Lactobacillaceae (12.059%, Phylum Firmicutes), and Barnesiellaceae (9.566%, Phylum Bacteroidetes). However, it is noteworthy that Peptostreptococcaceae is only present in four out of fourteen cloacal swab samples. These top three families found in cloacal data were present in all fourteen cecal content samples.
A PCoA (Fig. 3) was calculated to compare community member composition within the cecal content and cloacal swab methods. The cecal content samples cluster tightly together, whereas the cloacal swab samples show high variability while still encompassing the cecal content samples. This high variability is not surprising given the total number of families represented in the cloacal swabs. Overall, the ordination pattern of these paired samples shows broad-ranging differences between the two sampling approaches.
Finally, when comparing the cumulative distribution functions filtered for > 0% of the relative abundances between cecal and cloacal microbiotas, we found major differences in the structure of abundance. The KS (D = 0.11745, p-value = 4.692e-07) and Cucconi tests (C = 298.959, p-value = 0) were highly significant. The results of the Kolmogorov-Smirnoff and Cucconi tests on the cumulative distribution function further demonstrate that both the location and the scales of the cecal content and cloacal swab distributions of relative abundances are highly different.
3.3 Richness and diversity differences between cecal and cloacal samples
We compared microbial species richness and diversity of the cecal content and cloacal swabs using the Chao1 and Inverse Simpson estimators (Fig. 4). Both the Wilcoxon Signed-Rank test (W = 4, p-value = 0.000845) and paired t-test (t = 4.042, p-value = 0.001397) showed highly significant differences in the Chao1 index between the cecal content and cloacal swabs, with the highest richness observed in the cecal samples. A higher Chao1 value indicates a higher number of low abundance taxa, e.g., singletons [51, 52]. The higher value in cecal samples, suggests that more rare taxa were captured in cecal samples. However, the Inverse Simpson species diversity estimator was not different between cloacal and cecal samples based on a Wilcoxon Signed-Rank test (W = 36, p-value = 0.3258) and a paired t-test (t = 0.89623, p-value = 0.3864). Similar to the Chao1 findings, the cecal content had higher microbial diversity, compared to the cloacal swabs. As the Inverse Simpson index estimates the richness weighted by the proportional abundance of taxa present within samples, the non-significance suggests that the two types did not differ in their internal weighted abundances.
To assess whether cecal or cloacal swabs captured differences between dietary treatments, we performed richness and diversity analyses, comparing the two diets. Figure 5A shows comparisons within cecal data, and Fig. 5B shows cloacal swab comparisons. The Wilcoxon Signed-Rank test on the Chao1 estimator returned non-significant p-values for the cecal content treatments (W = 10, p-value = 0.1119) and cloacal swab treatments (W = 12, p-value = 0.1898). The Chao1 values were higher for T1 than T2 in both the cecal content and cloacal swab methods. Overall, neither the cecal nor the cloacal swab microbiotas differed between the dietary treatments based on the Chao1 index.
On the other hand, the Inverse Simpson index was not different based on cecal samples (W = 24, p-value = 0.8981), but the cloacal swab data showed significant differences between diets (W = 6, p-value = 0.02897). In summary, the cloacal samples showed significant differences between dietary treatments despite the absence of such distinction when comparing cecal microbiotas.
3.4 Cloacal swabs do not reflect the community structure inferred from cecal samples
The PERMANOVA analysis showed significant difference (F. Model = 4.903, R2 = 0.15866, p-value = < 0.001) in the centroids and dispersion between cecal and cloacal microbiota. This difference of community structure was further supported with significant results from the HOMOVA (BValue = 3.07143, p-value < 0.001) and Weighted Unifrac (WScore = 0.664192, WSig = < 0.001) tests. On the whole, the results from the PERMANOVA, HOMOVA, and Weighted Unifrac tests all show that the microbiota communities inferred from cloacal swabs are different from cecal microbiota.
Next, we investigated whether differences in the diet treatments (T1 and T2) elicit differences in communities (β diversity) inferred using cecal versus cloacal samples. We found that neither sampling method detected differences in β diversity between the diets. The results from PERMANOVA showed that the cecal content was not different between dietary treatments (F. Model = 1.4155, R2 = 0.10551, p-value = 0.1455). The cloacal swabs also revealed the same information (F. Model = 1.065, R2 = 0.08148, p-value = 0.3523). HOMOVA and AMOVA both yielded the same information, suggesting that between dietary treatments, the differences in variances (BValue = 0.518171, p-value = 0.125) or the genetic diversity of communities (BValue = 0.024201, p-value = 0.492) were not different. In contrast, the Weighted Unifrac comparing the cecal microbiota between dietary treatments was significant (WScore = 0.584437, WSig < 0.001), as was the comparison of cloacal microbiotas between diets (WScore = 0.77311, WSig < 0.001). However, in this instance, it is difficult to determine how the high variability among cloacal samples influence the comparisons of β diversity between dietary treatments.
Further comparison of the dietary treatments within each method and between methods using the KS, Lepage, and Cucconi tests yielded no significant differences (Table 2). We found non-significant results using the same three tests for the comparison of dietary treatments within each sampling method. These results are expected as we found no major difference in the structure and abundances between the treatments using either cecal content and cloacal swabs.
Table 2
Statistical testing results for geometric mean distributions of dietary treatments within and between methods
Comparison | KS Test | Lepage Test | Cucconi Test |
Cecal Content T1 vs. Cecal Content T2 | D = 0.3348; p-value = 0.306 | L = 2.777; p-value = 0.2555 | C = 1.394; p-value = 0.25 |
Cloacal Swab T1 vs. Cloacal Swab T2 | D = 0.3375; p-value = 0.2152 | L = 3.02; p-value = 0.2215 | C = 1.584; p-value = 0.194 |
Table 2. Summary of the statistical test used in comparing the geometric mean distributions. Geometric means of relative abundance were calculated and filtered for relative abundances > 2%. Specific dietary comparisons are in the first column, with the statistical tests in the remaining columns. Each cell contains the respective test statistic and p-value for that comparison.