Characterization of captive and wild 13-lined ground squirrel cecal microbiotas using Illumina-based sequencing
Animal Microbiome volume 4, Article number: 1 (2022)
Hibernating animals experience extreme changes in diet that make them useful systems for understanding host-microbial symbioses. However, most of our current knowledge about the hibernator gut microbiota is derived from studies using captive animals. Given that there are substantial differences between captive and wild environments, conclusions drawn from studies with captive hibernators may not reflect the gut microbiota’s role in the physiology of wild animals. To address this, we used Illumina-based sequencing of the 16S rRNA gene to compare the bacterial cecal microbiotas of captive and wild 13-lined ground squirrels (TLGS) in the summer. As the first study to use Illumina-based technology to compare the microbiotas of an obligate rodent hibernator across the year, we also reported changes in captive TLGS microbiotas in summer, winter, and spring.
Wild TLGS microbiotas had greater richness and phylogenetic diversity with less variation in beta diversity when compared to captive microbiotas. Taxa identified as core operational taxonomic units (OTUs) and found to significantly contribute to differences in beta diversity were primarily in the families Lachnospiraceae and Ruminococcaceae. Captive TLGS microbiotas shared phyla and core OTUs across the year, but active season (summer and spring) microbiotas had different alpha and beta diversities than winter season microbiotas.
This is the first study to compare the microbiotas of captive and wild rodent hibernators. Our findings suggest that data from captive and wild ground squirrels should be interpreted separately due to their distinct microbiotas. Additionally, as the first study to compare seasonal microbiotas of obligate rodent hibernators using Illumina-based 16S rRNA sequencing, we reported changes in captive TLGS microbiotas that are consistent with previous work. Taken together, this study provides foundational information for improving the reproducibility and experimental design of future hibernation microbiota studies.
Host-microbe symbioses are dynamic, especially in animals that experience extreme shifts in diet as these changes dramatically alter substrate availability for their gut microbiota. An example includes hibernating mammals. Hibernation is an ecophysiological strategy to survive times of reduced resource availability and high-energy demand by altering both behavior and physiology. In small mammalian hibernators such as the 13-lined ground squirrel (TLGS; Ictidomys tridecemlineatus), the circannual hibernation cycle involves periods of summer hyperphagia when the host acquires adequate fat stores for energy usage during winter hibernation . During summer, the gut microbiota has access to both dietary and host-derived substrates as rich sources of energy. In winter, the host enters hibernation, fasts for several months, and relies entirely upon fat stores for energy. The winter hibernation season is characterized by cycling between periods of depressed and normal metabolism (torpor and interbout arousal, respectively). Torpor occurs when metabolism slows to < 4% of normal rates, causing body temperature (Tb) to plummet to just above ambient temperature (< 10 °C) [2, 3]. Torpor is interrupted by brief interbout arousals (IBAs) that can last 12–24 h [2,3,4]. During IBAs, metabolism increases to normal rates and Tb warms to ~ 36 °C [2,3,4]. Due to fasting during hibernation, the gut microbiota is forced to rely solely on host-derived substrates (e.g., mucins) as its source of energy. The hibernation season ends with the emergence aboveground in spring and the resumption of normal metabolic activity and feeding patterns. This natural cycle of extreme changes in diet and physiology makes mammalian hibernators like the TLGS useful systems for studying host-microbe symbioses.
The majority of our current knowledge of the hibernator gut microbiota comes primarily from studies of captive animals. Captive environments differ from natural, wild environments due to necessary changes in diet, habitat, rearing conditions, and exposure to environmental microbes ; as a result, the microbiotas of captive animals likely differ substantially from those of wild animals. While the use of captive animals allows for controlled experiments that have led to significant advances in our understanding of their biology, their use for gut microbiota studies reduces confidence in the conclusions that can be drawn about the microbiota’s role in the physiology of wild animals. Moreover, these limitations may restrict the ability to translate discoveries using captive host-microbial relationships to other systems, such as the use of captive mice as a model for human health [5,6,7,8]. Differences in diet are of particular concern as diet is a known major driver of microbiota composition and metabolism [9,10,11,12,13,14,15,16,17] and it is difficult to recapitulate a hibernator’s wild diet in a captive setting. Additionally, the development of an animal’s microbiota early in life can have long-lasting effects on microbiota composition and host physiology [13, 18, 19], and many hibernation studies use animals born in the lab rather than in the wild. Although there is a growing number of studies comparing the microbiotas of captive and wild animals [5,6,7,8, 20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55], very few have been conducted with hibernating species. One study that examined three species of hibernating bats found that the number of operational taxonomic units (OTUs) is higher in captive than in wild bats . Studies in other animals have reported that differences between captive and wild microbiotas are species-specific [5, 30, 36, 53]. Therefore, it is important to understand how captivity alters the microbiota of hibernating animals like the TLGS.
Most rodent hibernation microbiota studies to date have used either 454 pyrosequencing [56,57,58,59] or clone libraries . However, advances in sequencing technology have resulted in the vast majority of microbiota studies using Illumina-based sequencing technology. Because different sequencing technologies may introduce biases, it is important to compare results among sequencing methods to understand how to interpret past data in the context of newer studies. Previous comparisons of 454 pyrosequencing and Illumina sequencing data demonstrate that the two methods result in similar beta diversity results [61,62,63]. This is also the case for alpha diversity for samples with high diversity but not for those with low diversity [62, 63]. Results from these two sequencing platforms also produce differences in microbiota taxonomic composition, with some studies reporting differences at the genus level  while others reveal differences at the phylum level . Therefore, it is important to consider how sequencing results from different platforms may differ in order to improve data interpretation.
Here, we used Illumina 16S rRNA gene sequencing to compare the bacterial cecal microbiotas of captive and wild TLGS in the summer. We hypothesize that wild TLGS microbiotas are more diverse than captive TLGS microbiotas, likely because the wild diet is more diverse than the captive diet. As this study is, to our knowledge, the first analysis of Illumina-based sequencing from an obligate rodent hibernator across active and hibernation seasons, we also reported changes in captive TLGS microbiotas across the year (summer, winter, and spring).
We sequenced a total of 59 samples, which consisted of 30 cecal content and 29 cecal mucosa. Bacterial amplicon sequencing generated a total of 3,252,557 raw sequences with an average of 54,209 ± 16,411 sequences per sample (mean ± SE; range 534–972,021). Sequence clean-up in mothur resulted in a total of 2,211,103 sequences for an average of 36,852 ± 13,449 sequences per sample (range 291–796,002). After normalization, 57 samples remained (29 content, 28 mucosa) with a range of reads = 10,453–11,159. All samples had Good’s coverage ≥ 95%, which indicated sufficient sequence coverage .
Cecal content and mucosa microbiotas are similar
To determine whether cecal content and mucosa microbiotas differ, we used paired tests to examine alpha and beta diversity within each of the five experimental groups. There were no differences in the number of observed OTUs, Faith’s phylogenetic diversity, phylogenetic evenness, or Shannon’s diversity (all paired t-test adj P ≥ 0.889, Additional file 1: Fig. S1). Similarly, we found no differences in the variances of weighted UniFrac, unweighted UniFrac, and Bray–Curtis dissimilarity (all betadisper adj P ≥ 0.293; Additional file 1: Fig. S2). We then examined the centroids of the three beta diversity metrics by including both squirrel ID and sample type in PERMANOVA tests. Squirrel ID accounted for the paired nature of the data and was significant for all three metrics in the captive summer group (all PERMANOVA adj P = 0.002) and for Bray–Curtis dissimilarity in the IBA group (adj P = 0.006), but not in the remaining groups and metrics (Additional file 2: Table S6). Within these significant groups and metrics, squirrel ID explained more variation in beta diversity than did sample type (Additional file 2: Table S6; all squirrel ID R2 ≥ 0.224, all sample type R2 ≤ 0.056). After accounting for squirrel ID, there were no differences in sample type for the three metrics (Additional file 2: Table S6; all adj P ≥ 0.912). Because sample type did not significantly impact alpha and beta diversity, we present data from mucosa microbiotas here and included those from cecal content in the supplementary materials.
Captive and Wild microbiotas have different alpha and beta diversities
We examined alpha diversity (within-sample diversity) between summer Captive and Wild groups by using four metrics (Fig. 1). In mucosa microbiotas, the number of observed OTUs was higher in Wild compared to Captive (Fig. 1A; T-test adj P < 0.001). Faith’s phylogenetic diversity, which is positively correlated with richness, also revealed the same difference as seen in number of observed OTUs in Wild vs. Captive (Fig. 1B; adj P < 0.001). There were no differences in phylogenetic evenness (Fig. 1C; adj P = 0.303) or Shannon’s diversity (Fig. 1D; adj P = 0.614). The same results were seen in content microbiotas (Additional file 1: Fig. S3).
We examined beta diversity by using weighted UniFrac, unweighted UniFrac, and Bray–Curtis dissimilarity (Fig. 2). In mucosa samples, the variance of unweighted UniFrac differed between Captive and Wild groups (Fig. 2; betadisper adj P = 0.012), but no group differences were observed in weighted UniFrac or Bray–Curtis Dissimilarity (Fig. 2; both adj P = 0.241). However, the centroids in all three metrics were different (Fig. 2; all adj P = 0.004). Content microbiotas had different variances and centroids between Captive and Wild groups in all three metrics (Additional file 1: Fig. S4; all betadisper adj P = 0.006, all PERMANOVA adj P = 0.014).
Few OTUs drive beta diversity differences between Captive and Wild microbiotas
To identify specific taxa that contributed to the differences between Captive and Wild microbiotas, we first examined phyla relative abundances. There were 10 phyla with total relative abundances > 1% across all mucosa samples (Table 1). Dominant phyla included the Firmicutes and Bacteroidetes, with smaller contributions from the Cyanobacteria, Tenericutes, Verrucomicrobia, Proteobacteria, Epsilonbacteraeota, Actinobacteria, Elusimicrobia, and Kiritimatiellaeota. We did not detect the Epsilonbacteraeota or Kiritimatiellaeota (Additional file 2: Table S7) in content microbiotas. Mucosa and content microbiotas had similar phyla relative abundances between Captive and Wild microbiotas (all adj P > 0.050); however, the phyla Proteobacteria and Cyanobacteria displayed a trend in content microbiotas (both adj P = 0.055).
As we saw no differences at the phylum-level, we proceeded at the OTU-level. Analysis of similarity percentages (SIMPER) in mucosa microbiotas identified three OTUs that contributed to differences in Bray–Curtis dissimilarity (Table 2). They explained 4.716% (cumulative SIMPER) of the difference in Bray–Curtis dissimilarity. Two of these OTUs belonged to the family Lachnospiraceae and one to the genus Lactobacillus. Both Lachnospiraceae OTUs had higher relative abundances in Wild compared to Captive, whereas the opposite was true for the Lactobacillus OTU (Table 2).
Content microbiotas had four OTUs that explained 7.260% of the Bray–Curtis dissimilarity (Additional file 2: Table S8). Three of these OTUs belonged to the family Lachnospiraceae, one of which was further classified to the genus Lachnospiraceae NK4A136 group, and they were not detected in Captive microbiotas (Additional file 2: Table S8). The remaining OTU was classified to the genus Lactobacillus and had higher relative abundance in Captive than in Wild (Additional file 2: Table S8). Two OTUs (one Lachnospiraceae and one Lactobacillus) were significant in both content and mucosa microbiotas.
Wild microbiotas have more core OTUs than Captive microbiotas
The Wild mucosa microbiota had 323 core OTUs whereas the Captive mucosa microbiota had 43 (Additional file 2: Table S9). The majority of the Wild core OTUs were classified to the families Ruminococcaceae (100 OTUs), Lachnospiraceae (93 OTUS), Muribaculaceae (34 OTUs), unclassified Gastranaerophilales (17 OTUs), and Rikenellaceae (14 OTUs). The majority of the Captive core OTUs were classified to the families Ruminococcaceae (15 OTUs) and Lachnospiraceae (10 OTUs). There were 24 shared core OTUs across all mucosa samples (Additional file 2: Table S9). All but one of these were classified to the order Clostridiales and the most common family was the Ruminococcaceae with 10 OTUs. The sole non-Clostridia OTU was classified to the genus Parasutterella. Content microbiotas had similar numbers and taxonomic classifications of core OTUs (Additional file 2: Table S10).
Captive microbiotas across seasons have similar phyla and core OTUs
Across all mucosa samples in the four captive groups, we identified 10 phyla with total relative abundances > 1% (Table 3). Only the Firmicutes differed among groups, with higher relative abundances in Summer compared to Torpor or IBA (both Tukey’s HSD adj P < 0.001). Phyla with no differences among groups included the Bacteroidetes, Verrucomicrobia, Proteobacteria, Kiritimatiellaeota, Cyanobacteria, Tenericutes, Actinobacteria, Elusimicrobia, and Epsilonbacteraeota. Content microbiotas had the same phyla, but we did not detect the Epsilonbacteraeota (Additional file 2: Table S11). In these samples, both the Firmicutes and Bacteroidetes differed among groups. The Firmicutes had higher relative abundances in Summer compared to Torpor or IBA (Additional file 2: Table S11; Dunn’s test adj P = 0.0012 and 0.001, respectively); the opposite was true for the Bacteroidetes (Additional file 2: Table S11; Tukey’s HSD Summer and Torpor adj P = 0.006, Summer and IBA adj P = 0.017). No other phyla differed among groups.
In the captive mucosa microbiotas, we identified 43 core OTUs in Summer, 65 in Torpor, 60 in IBA, and 141 in Spring (Table 4). Core OTUs in each group were dominated by the families Lachnospiraceae and Ruminococcaceae (Additional file 2: Table S12). Torpor and Spring also had > 10 core OTUs classified to the family Muribaculaceae. There were 16 core OTUs found in every mucosa sample. The content microbiotas had similar results (Additional file 2: Table S13), and 10 core OTUs found in all content samples were also found in all mucosa samples.
Captive microbiotas in active and winter groups have different alpha and beta-diversities
Summer mucosa microbiotas had more observed OTUs (Fig. 3A) compared to Torpor or IBA (Dunn’s test adj P = 0.006 and 0.010, respectively), but did not differ from Spring (adj P = 0.450). Torpor was not different from IBA or Spring (adj P = 0.408 and 0.043, respectively), and IBA was not different from Spring (adj P = 0.039). Phylogenetic diversity (Fig. 3B) was higher in Summer than IBA (adj P = 0.018), whereas phylogenetic evenness (Fig. 3C) was lower in Summer than Torpor (Tukey’s HSD adj P = 0.022). Shannon’s diversity did not differ among groups (Fig. 3D, ANOVA adj P = 0.868). Content microbiota results (Additional file 1: Fig. S5) were the same as for the mucosa microbiota with respect to the number of observed OTUs and Shannon’s diversity index; however, phylogenetic diversity and phylogenetic evenness differed. Phylogenetic diversity only displayed a trend of differences among groups (ANOVA adj P = 0.069). Phylogenetic evenness was lower in Summer compared to Torpor or IBA (Tukey’s HSD both adj P < 0.001), and lower in Spring compared to Torpor (adj P = 0.042).
Beta diversity variance across all captive mucosa microbiotas (Fig. 4) did not differ in weighted and unweighted UniFrac and in Bray–Curtis dissimilarity (all betadisper P ≥ 0.822). However, the centroids of active groups (Summer and Spring) were different from those of winter groups (Torpor and IBA) in weighted UniFrac (all PERMANOVA adj P < 0.005) and in Bray–Curtis dissimilarity (all adj P < 0.027). The unweighted UniFrac centroids differed between Summer and winter groups (Torpor adj P = 0.008, IBA adj P = 0.003), and between Spring and IBA (adj P = 0.033). There were no differences within active and winter groups (both adj P for weighted UniFrac ≥ 0.143, unweighted UniFrac ≥ 0.635, and Bray–Curtis dissimilarity ≥ 0.094). Content microbiota beta diversity results (Additional file 1: Fig. S6) were similar to those of mucosa microbiotas, except that the unweighted UniFrac centroids of Torpor and Spring were different (adj P = 0.038). Lastly, we attempted to identify OTUs that drove these differences in beta diversity by using SIMPER but detected no significant OTUs (all adj P ≥ 0.075).
In this study, we sought to determine if the bacterial cecal microbiotas of TLGS were impacted by captive versus wild environments. We compared summer microbiotas between TLGS born in captivity that consumed a chow diet with wild-caught TLGS that consumed a natural, wild diet. We found that Wild microbiotas had greater richness and phylogenetic diversity and decreased variance in beta diversity compared to Captive microbiotas. Important taxa that contributed to differences in beta diversity and were core OTUs were primarily classified to the families Lachnospiraceae and Ruminococcaceae. We also described the microbiotas of captive TLGS across the year (summer, winter, and spring) using Illumina-based sequencing as most previous rodent hibernation microbiota studies to date have used either 454 pyrosequencing [56,57,58,59] or clone libraries .
Animals born and raised in the wild experience several environmental features that can affect their gut microbiota composition and therefore influence the host-microbe symbiosis. The composition and availability of food is arguably the most important environmental influence as these can affect early microbiota development and host-microbe relationships throughout life [12, 13, 16, 17]. Although commercial animal chow provides adequate macro- and micro-nutrients for generalized rodents, precise recapitulation of the wild diet is difficult, especially for omnivores like the TLGS. Foraging in the wild also provides opportunities for animals to ingest microbes that contribute to their gut microbiota. This variable is largely missing in captive diets, which can have significant effects on microbiota diversity and the host-microbe symbiosis .
Another potential variable is host age. Studying animals in captivity is useful because it provides precise age information that is often unavailable for wild animals. In our study, we were unable to assign exact ages to the Wild group; however, based on their body weights at capture (Additional file 2: Table S1), we predicted that they were likely adults of at least 13 months of age. Therefore, the Wild squirrels were all older than the Captive squirrels. Previous TLGS microbiota studies have found that there are few differences between the summer microbiotas of wild-caught mothers and the summer microbiotas of their captive pups that were born and raised in captivity with known ages [56, 57]. Therefore, we posit that age differences between our Wild and Captive TLGS would not significantly bias our results, but future studies should account for this variable.
Additional limitations of our study include the use of small sample sizes and relative abundance data. Due to our small sample sizes, we only have sufficient power to detect very large effect sizes that were often larger than our observed effect sizes (Additional file 2: Tables S4 and S5) [65, 66]. Although the power of our study is limited, our work still provides crucial information to help future studies calculate the number of samples needed to have sufficient statistical power. Future studies can reference our observed effect sizes to determine sample sizes for both wild and captive ground squirrels, whereas such data was not previously available. Additionally, our use of relative abundance data may be misleading as changes in the abundance of one taxa will influence the relative abundances of all other taxa in that sample. This is especially important in hibernating ground squirrels because absolute bacterial abundance dramatically decreases in hibernation compared to active seasons .
Overall, we found that differences in the alpha diversity of content and mucosa microbiotas between the Captive and Wild groups were attributed to the Wild group’s larger number of OTUs with low relative abundances. The Wild group had significantly more observed OTUs (Fig. 1A) and higher phylogenetic diversity (Fig. 1B) compared to Captive. These additional Wild OTUs were from numerous different taxa rather than from a single or closely related group of taxa as there was no difference in the average phylogenetic distance among all species pairs (Fig. 1C). Most of the Wild OTUs also had low relative abundances, demonstrated by the lack of difference in Shannon’s index (Fig. 1D).
Beta diversity of content and mucosa microbiotas between Wild and Captive groups were different with respect to OTU identities, relative abundances, and phylogenetic relatedness (Fig. 2). Wild content microbiotas also had less variation than Captive content microbiotas; in mucosa microbiotas, this was true only in unweighted UniFrac. Because diet is a known driver of the microbiota, we expected the Wild group to have higher variation compared to the Captive group, as reported in a previous study that examined three species of hibernating bats . However, we found the opposite to be true. This was consistent with another study that found captive rodent microbiotas had more variation than wild rodent microbiotas . One possible explanation is that the wild squirrels were captured from the same area on the same day, so their diets may have been very similar. Additionally, the small sample size of the Wild group may underestimate beta diversity variation.
Although Captive and Wild microbiotas had different alpha and beta diversities, we detected the same phyla in both groups, with no significant differences in relative abundances (Table 1); however, we did identify trends in the Proteobacteria and Cyanobacteria in content microbiotas. The Proteobacteria displayed a trend of higher relative abundance in Captive than in Wild microbiotas, which contrasted with results in other wild rodent studies where the Proteobacteria have higher relative abundances in wild compared to captive animals [6, 42, 67]. The top Proteobacteria OTUs in Captive microbiotas included members of the family Desulfovibrionaceae, such as the genera Desulfovibrio and Parasutterella. Desulfovibrio is a sulfate-reducer known to participate in mucin glycan degradation  and increases in relative abundance during fasting  and hibernation [56, 59]. Parasutterella isolates from rodent and human guts are asaccharolytic and predicted to rely on amino acids to survive ; they also induce changes in bile acid metabolism and aromatic amino acid catabolism when introduced in mice . The Cyanobacteria had higher relative abundance in Captive than in Wild microbiotas. All TLGS Cyanobacteria OTUs were classified to the class Melainabacteria and the order Gastranaerophilales. Melainabacteria, a proposed sister phylum to Cyanobacteria, is reported to not be autotrophic [70,71,72], has been detected in diverse environments ranging from the human gut to groundwater , and is predicted to produce formate and ethanol in the gut .
We identified the families Lachnospiraceae and Ruminococcaceae as the predominant classifications of core OTUs in both Captive and Wild microbiotas (Additional file 2: Tables S5 and S6) and of OTUs that explained Bray–Curtis dissimilarity (Table 2 and Additional file 2: Table S4). Members of both families are known to degrade diverse plant polysaccharides [73, 74]. There were more core OTUs in Wild compared to Captive microbiotas, which is not surprising given that Wild microbiotas had more observed OTUs; however, this could also be due to its small sample size (n = 4). Just under half of the Captive core OTUs were also Wild core OTUs, demonstrating that there is reasonable overlap between the two groups despite their different alpha and beta diversities. Bray–Curtis dissimilarity between Captive and Wild microbiotas was explained by several Lachnospiraceae OTUs that had higher relative abundances in Wild than in Captive microbiotas, which is consistent with what has been observed in wild and captive mice [31, 39]. These taxa likely help the host metabolize plant material, as members of the Lachnospiraceae are known for degrading diverse polysaccharides and are positively correlated with acetate, butyrate, and propionate in TLGS .
It is important to note that most rodent hibernation microbiota studies have used either 454 pyrosequencing [56,57,58,59] or clone libraries . However, advances in sequencing technology have resulted in the vast majority of microbiota studies using Illumina-based sequencing. Because different sequencing technologies are known to introduce various biases [61,62,63], we considered if seasonal changes in the microbiotas of captive TLGS across the year (summer, winter, and spring) were consistent between Illumina sequencing and pyrosequencing. Overall, our findings from Illumina-based sequencing were generally consistent with those based on 454 pyrosequencing. We identified more phyla than in previous studies, but only found changes in the relative sequence abundances of the dominant phyla, and not in those with lower relative abundances [56, 57]. This suggests that trends in taxa with high relative abundances are conserved regardless of sequencing technology, but there may be differences in taxa with low relative abundances. These differences may also be due to different primers used in one of the previous studies . We also confirmed previous findings that a small number of core OTUs persist in the TLGS cecal microbiota throughout the year despite seasonal shifts in the host’s diet . Finally, our results with respect to alpha and beta diversity in winter (Torpor and IBA) and active groups (Summer and Spring) were similar, although we note that there were some minor differences in the gut microbiotas of active group squirrels. Taken together, these data confirm previous findings on the TLGS gut microbiota and address concerns that may arise due to biases inherent to the two sequencing technologies.
This is the first study to compare the microbiotas of captive and wild rodent hibernators. We demonstrated that Wild TLGS microbiotas have increased richness and phylogenetic diversity, and decreased variation in beta diversity compared to Captive TLGS microbiotas. In both Captive and Wild microbiotas, important taxa that are core OTUs or significantly contribute to beta diversity are predominantly in the families Lachnospiraceae and Ruminococcaceae. Given the significant differences between wild and captive TLGS microbiotas, it is important that future research determines how long it takes for wild microbiotas to change in captivity and what those changes entail. We also reported seasonal differences in captive TLGS microbiotas throughout the annual cycle by using Illumina-based 16S rRNA gene sequencing. Results about alpha diversity, beta diversity, and taxa with high relative abundances were consistent with past conclusions based on 454 pyrosequencing, but some differences emerged when comparing Summer and Spring seasons and in the analysis of low abundance taxa. Because next-generation sequencing continues to advance quickly, continued methodological comparisons should be conducted to evaluate any biases introduced by new methods. Taken together, our results help improve reproducibility and experimental design of future hibernation microbiota studies.
All procedures were approved by the University of Wisconsin-Madison Institutional Animal Care and Use Committee under protocols V001229 and V005481. Two groups of healthy TLGS were used in our study: a captive and a wild group (Fig. 5). The captive group consisted of 26 animals (17 females and 9 males) that were born in captivity to seven pregnant, wild-caught females in Madison, WI. The pregnant females were housed individually at 22 °C with a 12:12 h light–dark cycle, provided water and rat chow (Harlan Teklad no. 7001, Indianapolis, IN, USA) ad libitum, and supplemented with apples, strawberries, and sunflower seeds weekly. Pups were born in May 2016 and remained with their mothers for five weeks before separation into new cages that housed two pups per cage. After two weeks, pups were transferred into individual cages. Water and chow were provided ad libitum for two weeks. Then chow was restricted to 12 g/day to prevent excessive weight gain and sunflower seeds (~ 1 g) were provided weekly. The wild squirrel group consisted of four adult TLGS (two males and two females) captured in July 2017 in Madison, WI.
The 26 captive TLGS were randomly assigned to one of four groups (Fig. 5): Summer (n = 8), Torpor (n = 6), IBA (n = 8), or Spring (n = 4). Captive Summer squirrels were sampled in August 2016. The remaining captive animals were transferred to a 4 °C room in September/October 2016 for hibernation. The room was held in constant darkness except for a daily, brief period of dim light (~ 5 min) to check activity states using the sawdust method . This involved placing very small pieces of paper on top of a torpid squirrel and examining if the paper had been disturbed, which would indicate that a squirrel had entered an IBA and then reentered torpor since the previous check. Once squirrels began using torpor, food and water were removed. Winter squirrels were sampled in January–February (after ~ 3.5–4.5 months in the cold room). Torpor squirrels were sampled during torpor (Tb < 10 °C). For IBA squirrels, torpid animals were brought to a lit 22 °C room for three hours to induce an arousal and sampled when Tb > 34 °C. Spring squirrels were removed from the cold room in January 2017 and transferred to a warm room with food and water ad libitum for two weeks before sampling in February 2017. Samples were collected from wild summer squirrels on the same day they were captured in July 2017. Metadata for all squirrels are listed in Additional file 2: Table S1.
Four anesthesia/euthanasia methods were used due to equipment availability and sample collection requirements for a separate project. Captive Summer and Spring groups were euthanized via exposure to isoflurane for 5–15 min followed by decapitation. Torpor squirrels were euthanized by decapitation or cervical dislocation. IBA squirrels were anesthetized via exposure to isoflurane or CO2 for 5–15 min followed by decapitation (Additional file 2: Table S1). Wild summer squirrels were euthanized with pentobarbital followed by decapitation. To determine whether anesthesia method had a significant impact on the gut microbiota, IBA squirrels that received isoflurane or CO2 were compared as this was the only experimental group where multiple anesthesia methods were used (Additional file 2: Table S2). We did not find evidence that the microbiotas of isoflurane- or CO2-administered animals were different. The methods [76,103,78] and results for this analysis are shown in Additional file 2: Table S2.
Tb was measured immediately by inserting a clean thermal probe into the body cavity. Cecal contents were collected by removing intact ceca and gently scraping the content into a sterile tube. To remove any remaining content, cecal tissue was rinsed with sterile phosphate-buffered saline (PBS) and any remaining liquid was gently scraped off. Cecal mucosae were collected using a razor blade to gently scrape off the mucosa. To assess the captive diet microbiota, three chow samples were collected. A representative wild diet was not collected as TLGS are omnivorous and consume a diverse diet that includes various plants, insects, and small birds . All samples were placed on dry ice and stored at − 80 °C. We collected a total of 62 samples: 30 cecal content samples (26 captive, 4 wild), 29 cecal mucosa samples (25 captive, 4 wild), and three chow samples.
Total genomic DNA was extracted from each sample using a phenol:chloroform extraction protocol  with the following modification: all aqueous phase washes used 25:24:1 phenol:chloroform:isoamyl alcohol instead of phenol:chloroform for a total of four washes. Four controls were processed with the cecal samples. Two sample collection controls were aliquots from the sterile PBS used to rinse cecal tissue. These were uncovered and exposed to air for the same amount of time needed to collect the cecal samples. Two extraction method controls contained sterile water. DNA was quantified with the Qubit Fluorometer (Invitrogen, San Diego, CA, USA) and stored at − 80 °C.
DNA amplification and sequencing
We used universal bacterial primers flanking the V4 region of the 16S rRNA gene (515F: GTGCCAGCMGCCGCGGTAA, 806R: GGACTACHVGGGTWTCTAAT) [81, 82]. The primers also contained adapters (forward: AATGATACGGCGACCACCGAGATCTACAC, reverse: CAAGCAGAAGACGGCATACGAGAT) compatible with Illumina sequencing technology (Illumina, San Diego, CA, USA) and unique barcodes for multiplexing (forward: eight unique eight bp barcodes, reverse: 12 unique eight bp barcodes). Each reaction contained 50 ng DNA, 0.4 μM forward primer, 0.4 μM reverse primer, 12.5 μL 2X HotStart ReadyMix (KAPA Biosystems, Wilmington, MA, USA), and water to a final volume of 25 μL. Polymerase chain reaction (PCR) was performed using a Bio-Rad S1000 thermocycler (Bio-Rad Laboratories, Hercules, CA, USA). Cycling conditions began with initial denaturation at 95 °C for 3 min, followed by 25 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and a final extension at 72 °C for 5 min. Four controls consisting of sterile water were processed alongside sample DNA to ensure there was no contamination during PCR. PCR products were purified using gel extraction from a 1.0% low-melt agarose gel (National Diagnostics, Atlanta, GA, USA) with a ZR-96 Zymoclean DNA Recovery Kit (Zymo Research, Irvine, CA, USA) and DNA was quantified using a Qubit Fluorometer. 60 samples were successfully amplified, and two samples failed to amplify (both chow samples). Samples were equimolarly pooled with 10% PhiX control DNA and sequenced on an Illumina MiSeq using a MiSeq 2 × 250 v2 kit.
Microbiota sequence clean-up
Sequences were demultiplexed on the Illumina MiSeq and processed using mothur v1.44.3  and SILVA release 132 [84,82,83,87]. Fastq files were submitted to NCBI’s Short Read Archive and are publicly available under accession number PRJNA742778. Mothur logfiles and output files are found at https://github.com/ednachiang/CaptiveWild/tree/master/mothur_output.
We followed the mothur standard operating procedure accessed on March 19, 2021 with the following modifications: we used screen.seqs maxlength = 300 and chimera.uchime. Briefly, the procedure involved combining paired reads, removing poor quality reads, aligning reads to SILVA release 132 [84,82,83,87], and removing undesirable reads (chimeras, Archaea, Eukaryota, chloroplasts, mitochondria, and unknown). Sequences were grouped into 97% operational taxonomic units (OTUs). Samples were normalized to 11,013 sequences using a modified total sum scaling method that optimized sample coverage. First, the proportion of each OTU relative to total reads was calculated separately for each sample. To ensure that all samples had sufficient coverage, a read cutoff was selected based on the sample with the lowest number of sequences that resulted in all normalized samples having Good’s coverage ≥ 95% . The OTU proportions were then multiplied by the sequence cutoff and rounded to the nearest integer. Three samples were removed due to insufficient sequencing (Good’s coverage < 95%): one chow, one Torpor mucosa, and one summer Wild content sample. All eight control samples were also discarded due to low number of sequences (mean = 32; range 3–117). To create a phylogenetic tree, a representative sequence from each OTU was first selected using get.oturep and the most abundant sequence. Representative sequences were renamed to match their assigned OTU and a phylip-formatted distance matrix was calculated. This distance matrix was used to calculate a phylogenetic tree using clearcut and the tree was classified using classify.tree.
Statistical analysis in R
Outputs from mothur were imported into R version 3.4.3  using the phyloseq  and phytools  packages. All visualizations were created using the ggplot2 , grid , and gridExtra  packages. The following packages were also used for statistical analyses: ape , dplyr , dunn.test , effsize , equivalence , picante , pwr , stats , tidyr , and vegan .
We used different subsets of samples to perform three analyses. To compare cecal content and mucosa microbiotas, only squirrels with both sample types were included due to the use of paired tests, and samples within each experimental group were analyzed separately. To compare Captive and Wild microbiotas, only summer Captive and summer Wild squirrels were considered, and content and mucosa samples were analyzed separately. Lastly, to compare captive microbiotas across the year (summer, winter, and spring), we analyzed captive groups (Summer, Torpor, IBA, and Spring) separately for content and mucosa samples. In all three comparisons, OTUs found in less than three samples were removed. Four outlier samples were also removed: content and mucosa samples from one IBA and one Spring individual (Additional file 2: Table S3). These were identified as samples with: (1) alpha diversity that was more than 1.5 times the interquartile range above the upper quartile or below the lower quartile for its experimental group, and (2) large beta diversity distances that, in an ordination, clearly clustered away from all other samples in the same experimental group. Additional information about the four outlier samples is found in Additional file 2: Table S3. All code to replicate statistical analyses and generate figures is available at https://github.com/ednachiang/CaptiveWild.
Alpha diversity analysis
We evaluated alpha diversity (within-sample diversity) using the number of observed OTUs, Faith’s phylogenetic diversity, phylogenetic evenness, and Shannon’s weighted diversity index . The number of observed OTUs was calculated using the phyloseq package . Faith’s phylogenetic diversity was computed by taking the sum of all branches in a tree using the picante package (pd) . Because phylogenetic diversity is positively correlated with richness, we included a metric unbiased by richness: mean pairwise distance (MPD), which measures phylogenetic evenness . MPD represents the average phylogenetic distance between all species pairs. Positive values indicate more phylogenetic evenness while negative values indicate less phylogenetic evenness. MPD was calculated by comparing phylogenetic diversity to a null model over 999 iterations (cophenetic, stats package ; ses.mpd, picante package ). Lastly, Shannon’s weighted diversity index was calculated in mothur . We tested the normality of each alpha diversity metric using quantile–quantile plots and the Shapiro test (qqnorm, qqline, shapiro.test, stats package ) before selecting a statistical test.
To compare content and mucosa microbiotas, paired t-tests were used as all metrics had normal distributions. Similarly, t-tests were used to compare Captive and Wild microbiotas. To compare captive microbiotas across seasons, content samples were analyzed using analysis of variance (ANOVA) (aov, stats package ) followed by pairwise Tukey’s Honest Significant Difference tests (TukeyHSD, stats package ). Mucosa MPD and Shannon’s diversity index were also analyzed this way. Mucosa number of observed OTUs and phylogenetic diversity were not normally distributed and were therefore analyzed using Kruskal–Wallis (kruskal.test, stats package ) followed by Dunn’s Test (dunn.test, dunn.test package ). All p-values were adjusted for false discovery rate using the Benjamini–Hochberg procedure (p.adjust, stats package ). Results were considered significant if adj P < 0.05 except for Dunn’s test where adj P ≤ 0.025 was significant.
Beta diversity analysis
Beta diversity (between-sample diversity) was examined using weighted UniFrac, unweighted UniFrac, and Bray–Curtis dissimilarity. All three metrics were calculated using the phyloseq package (distance)  and visualized using principal coordinate analysis (PCoA) ordinations (ordinate, phyloseq package ). To test whether there were significant differences in variance/dispersion, homogeneity of groups dispersions tests (betadisper, vegan package ) were used. Next, we tested if beta diversity centroids differed among experimental groups using permutational multivariate analysis of variance  (PERMANOVA; adonis, vegan package ). All p-values were corrected for false discovery rate using the Benjamini–Hochberg procedure.
To identify which OTUs significantly contributed to differences in Bray–Curtis dissimilarity, the similarity percentages test (SIMPER; simper, vegan package ) was used along with the Kruskal–Wallis test with false discovery rate correction using the Benjamini–Hochberg procedure. Only OTUs that accounted for ≥ 1% of the differences in beta diversity with adj P < 0.05 were considered.
Phylum-level relative abundances were calculated using the phyloseq  and dplyr  packages and phyla with total relative abundances < 1% were removed. Prior to selecting a statistical test, we evaluated normality using quantile–quantile plots and Shapiro tests. For the Captive and Wild comparison, phyla with normal distributions were analyzed using t-tests, whereas those with non-normal distributions were analyzed using Wilcoxon Rank Sum tests (wilcox.test, stats package ). To compare captive microbiotas across seasons, phyla with normal distributions were analyzed using ANOVA followed by Tukey’s HSD and those with non-normal distributions were analyzed using Kruskal–Wallis followed by Dunn’s test. All p-values were corrected for false discovery rate using the Benjamini–Hochberg procedure and results were considered significant if adj P < 0.05, except for Dunn’s test where adj P ≤ 0.025 was significant.
Identification of core OTUs
Core OTUs were defined as OTUs that were present in every sample within a group and were identified using the phyloseq package . Core OTUs shared among groups were identified by taking the intersect of the core OTUs in two or more groups.
Calculation of effect size and power
Effect size and power were calculated for alpha diversity comparisons between paired content and mucosa samples (Additional file 2: Table S4), and between summer Captive and Wild samples (Additional file 2: Table S5). They were also calculated for comparisons of phyla relative abundances between summer wild and captive samples (Additional file 2: Table S5). Effect size (Cohen’s d) was calculated with cohen.d (effsize package ), and power was calculated with pwr.t2n.test (pwr package ).
Availability of data and material
All FAST Q files were submitted to the NCBI’s Short Read Archive and are publicly available under bioproject number PRJNA742778 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA742778). Mothur logfiles and output files, as well as all R code used to perform statistical analyses and generate figures, are publicly available at https://github.com/ednachiang/CaptiveWild.
Analysis of variance
Mean pairwise distance
Operational taxonomic unit
Principal coordinate analysis
Polymerase chain reaction
Permutational multivariate analysis of variance
Similarity percentages test
- Tb :
13-Lined ground squirrel
Carey HV, Assadi-Porter FM. The hibernator microbiome: host-bacterial interactions in an extreme nutritional symbiosis. Annu Rev Nutr. 2017;37:477–500.
Johnson GE. Hibernation of the thirteen-lined ground squirrel, Citellus tridecemlineatus (Mitchell). J Exp Zool. 1928;50(1):15–30.
Carey HV, Andrews MT, Martin SL. Mammalian hibernation: cellular and molecular responses to depressed metabolism and low temperature. Physiol Rev. 2003;83:1153–81.
Johnson GE. Hibernation of the thirteen-lined ground squirrel, Citellus tridecemlineatus (Mitchill). III. The rise in respiration, heart beat and temperature in waking from hibernation. Biol Bull. 1929;57(2):107–29.
McKenzie VJ, Song SJ, Delsuc F, Prest TL, Oliverio AM, Korpita TM, Alexiev A, Amato KR, Metcalf JL, Kowalewski M, et al. The effects of captivity on the mammalian gut microbiome. Integr Comp Biol. 2017;57(4):690–704.
Rosshart SP, Vassallo BG, Angeletti D, Hutchinson DS, Morgan AP, Takeda K, Hickman HD, McCulloch JA, Badger JH, Ajami NJ, et al. Wild mouse gut microbiota promotes host fitness and improves disease resistance. Cell. 2017;171(5):1015–28.
Rosshart SP, Herz J, Vassallo BG, Hunter A, Wall MK, Badger JH, McCulloch JA, Anastasakis DG, Sarshad AA, Leonardi I, et al. Laboratory mice born to wild mice have natural microbiota and model human immune responses. Science. 2019;365(6452):461.
Yeung F, Chen YH, Lin JD, Leung JM, McCauley C, Devlin JC, Hansen C, Cronkite A, Stephens Z, Drake-Dunn C, et al. Altered immunity of laboratory mice in the natural environment is associated with fungal colonization. Cell Host Microbe. 2020;27(5):809–22.
Herter CA, Kendall AI. The influence of dietary alterations on the types of intestinal flora. J Biol Chem. 1910;7:203–36.
Cannon PR. The effects of diet on the intestinal flora. J Infect Dis. 1921;29(4):369–85.
Louis P, Scott KP, Duncan SH, Flint HJ. Understanding the effects of diet on bacterial metabolism in the large intestine. J Appl Microbiol. 2007;102(5):1197–208.
Turnbaugh PJ, Ridaura VK, Faith JJ, Rey FE, Knight R, Gordon JI. The effect of diet on the human gut microbiome: a metagenomic analysis in humanized gnotobiotic mice. Sci Transl Med. 2009;1(6):614.
Buddington RK, Sangild PT. Companion animals symposium: development of the mammalian gastrointestinal tract, the resident microbiota, and the role of diet in early life. J Anim Sci. 2011;89(5):1506–19.
David LA, Maurice CF, Carmody RN, Gootenberg DB, Button JE, Wolfe BE, Ling AV, Devlin AS, Varma Y, Fischbach MA, et al. Diet rapidly and reproducibly alters the human gut microbiome. Nature. 2014;505(7484):559–63.
Voreades N, Kozil A, Weir TL. Diet and the development of the human intestinal microbiome. Front Microbiol. 2014;5:494.
Kalmokoff M, Franklin J, Petronella N, Green J, Brooks SP. Phylum level change in the cecal and fecal gut communities of rats fed diets containing different fermentable substrates supports a role for nitrogen as a factor contributing to community structure. Nutrients. 2015;7(5):3279–99.
Krautkramer KA, Kreznar JH, Romano KA, Vivas EI, Barrett-Wilt GA, Rabaglia ME, Keller MP, Attie AD, Rey FE, Denu JM. Diet-microbiota interactions mediate global epigenetic programming in multiple host tissues. Mol Cell. 2016;64(5):982–92.
Backhed F, Roswall J, Peng Y, Feng Q, Jia H, Kovatcheva-Datchary P, Li Y, Xia Y, Xie H, Zhong H, et al. Dynamics and stabilization of the human gut microbiome during the first year of life. Cell Host Microbe. 2015;17(5):690–703.
Martinez I, Maldonado-Gomez MX, Gomes-Neto JC, Kittana H, Ding H, Schmaltz R, Joglekar P, Cardona RJ, Marsteller NL, Kembel SW, et al. Experimental evaluation of the importance of colonization history in early-life gut microbiota assembly. Elife. 2018;7:e36521.
Uenishi G, Fujita S, Ohashi G, Kato A, Yamauchi S, Matsuzawa T, Ushida K. Molecular analyses of the intestinal microbiota of chimpanzees in the wild and in captivity. Am J Primatol. 2007;69(4):367–76.
Scupham AJ, Patton TG, Bent E, Bayles DO. Comparison of the Cecal microbiota of domestic and wild Turkeys. Microb Ecol. 2008;56:322–31.
Villers LM, Jang SS, Lent CL, Lewin-Koh SC, Norosoarinaivo JA. Survey and comparison of major intestinal flora in captive and wild ring-tailed lemur (Lemur catta) populations. Am J Primatol. 2008;70(2):175–84.
Xenoulis PG, Gray PL, Brightsmith D, Palculict B, Hoppes S, Steiner JM, Tizard I, Suchodolski JS. Molecular characterization of the cloacal microbiota of wild and captive parrots. Vet Microbiol. 2010;146(3–4):320–5.
Dhanasiri AKS, Brunvold L, Brinchmann MF, Korsnes K, Bergh Ø, Kiron V. Changes in the intestinal microbiota of wild Atlantic cod Gadus morhua L. upon captive rearing. Microb Ecol. 2011;61:20–30.
Nakamura N, Amato KR, Garber P, Estrada A, Mackie RI, Gaskins HR. Analysis of the hydrogenotrophic microbiota of wild and captive black howler monkeys (Alouatta pigra) in palenque national park, Mexico. Am J Primatol. 2011;73(9):909–19.
Wienemann T, Schmitt-Wagner D, Meuser K, Segelbacher G, Schink B, Brune A, Berthold P. The bacterial microbiota in the ceca of Capercaillie (Tetrao urogallus) differs between wild and captive birds. Syst Appl Microbiol. 2011;34(7):542–51.
Eigeland K. Bacterial community structure in the hindgut of wild and captive dugongs (Dugong dugon). Aquat Mamm. 2012;38(4):402–11.
Nelson TM, Rogers TL, Carlini AR, Brown MV. Diet and phylogeny shape the gut microbiota of Antarctic seals: a comparison of wild and captive animals. Environ Microbiol. 2013;15(4):1132–45.
Kohl KD, Dearing MD. Wild-caught rodents retain a majority of their natural gut microbiota upon entrance into captivity. Environ Microbiol Rep. 2014;6(2):191–5.
Kohl KD, Skopec MM, Dearing MD. Captivity results in disparate loss of gut microbial diversity in closely related hosts. Conserv Physiol. 2014;2(1):cou009.
Kreisinger J, Cizkova D, Vohanka J, Pialek J. Gastrointestinal microbiota of wild and inbred individuals of two house mouse subspecies assessed using high-throughput parallel pyrosequencing. Mol Ecol. 2014;23(20):5048–60.
Clayton JB, Vangay P, Huang H, Ward T, Hillmann BM, Al-Ghalith GA, Travis DA, Long HT, Tuan BV, Minh VV, et al. Captivity humanizes the primate microbiome. PNAS. 2016;113(37):10376–81.
Delport TC, Power ML, Harcourt RG, Webster KN, Tetu SG. Colony location and captivity influence the gut microbial community composition of the Australian Sea Lion (Neophoca cinerea). Appl Environ Microbiol. 2016;82(12):3440–9.
Xie Y, Xia P, Wang H, Yu H, Giesy JP, Zhang Y, Mora MA, Zhang X. Effects of captivity and artificial breeding on microbiota in feces of the red-crowned crane (Grus japonensis). Sci Rep. 2016;6:33350.
Borbon-Garcia A, Reyes A, Vives-Florez M, Caballero S. Captivity shapes the gut microbiota of andean bears: insights into health surveillance. Front Microbiol. 2017;8:1316.
Kohl KD, Brun A, Magallanes M, Brinkerhoff J, Laspiur A, Acosta JC, Caviedes-Vidal E, Bordenstein SR. Gut microbial ecology of lizards: insights into diversity in the wild, effects of captivity, variation across gut regions and transmission. Mol Ecol. 2017;26(4):1175–89.
Li Y, Hu X, Yang S, Zhou J, Zhang T, Qi L, Sun X, Fan M, Xu S, Cha M, et al. Comparative analysis of the gut microbiota composition between captive and wild forest musk deer. Front Microbiol. 2017;8:1705.
Allan N, Knotts TA, Pesapane R, Ramsey JJ, Castle S, Clifford D, Foley J. Conservation implications of shifting gut microbiomes in captive-reared endangered voles intended for reintroduction into the wild. Microorganisms. 2018;6(3):94.
Leung JM, Budischak SA, Chung The H, Hansen C, Bowcutt R, Neill R, Shellman M, Loke P, Graham AL. Rapid environmental effects on gut nematode susceptibility in rewilded mice. PLoS Biol. 2018;16(3):e2004108.
Chong R, Grueber CE, Fox S, Wise P, Barrs VR, Hogg CJ, Belov K. Looking like the locals—gut microbiome changes post-release in an endangered species. Anim Microbiome. 2019;1(1):8.
Gibson KM, Nguyen BN, Neumann LM, Miller M, Buss P, Daniels S, Ahn MJ, Crandall KA, Pukazhenthi B. Gut microbiome differences between wild and captive black rhinoceros—implications for rhino health. Sci Rep. 2019;9(1):7570.
Schmidt E, Mykytczuk N, Schulte-Hostedde AI. Effects of the captive and wild environment on diversity of the gut microbiome of deer mice (Peromyscus maniculatus). ISME J. 2019;13(5):1293–305.
Shinohara A, Nohara M, Kondo Y, Jogahara T, Nagura-Kato GA, Izawa M, Koshimoto C. Comparison of the gut microbiotas of laboratory and wild Asian house shrews (Suncus murinus) based on cloned 16S rRNA sequences. Exp Anim. 2019;68(4):531–9.
Smith T. The effect of gut microbiota on host fitness of animals released from captivity. Groningen: University of Groningen; 2019.
Xiao Y, Xiao G, Liu H, Zhao X, Sun C, Tan X, Sun K, Liu S, Feng J. Captivity causes taxonomic and functional convergence of gut microbial communities in bats. PeerJ. 2019;7:e6844.
Edenborough KM, Mu A, Muhldorfer K, Lechner J, Lander A, Bokelmann M, Couacy-Hymann E, Radonic A, Kurth A. Microbiomes in the insectivorous bat species Mops condylurus rapidly converge in captivity. PLoS ONE. 2020;15(3):e0223629.
Lin JD, Devlin JC, Yeung F, McCauley C, Leung JM, Chen YH, Cronkite A, Hansen C, Drake-Dunn C, Ruggles KV, et al. Rewilding Nod2 and Atg16l1 mutant mice uncovers genetic and environmental contributions to microbial responses and immune cell composition. Cell Host Microbe. 2020;27(5):830–40.
Martinez-Mota R, Kohl KD, Orr TJ, Denise Dearing M. Natural diets promote retention of the native gut microbiota in captive rodents. ISME J. 2020;14(1):67–78.
Ni Q, He X, Zeng B, Meng X, Xu H, Li Y, Yang M, Li D, Yao Y, Zhang M, et al. Variation in gut microbiota of captive bengal slow lorises. Curr Microbiol. 2020;77(10):2623–32.
Ning Y, Qi J, Dobbins MT, Liang X, Wang J, Chen S, Ma J, Jiang G. Comparative analysis of microbial community structure and function in the gut of wild and captive Amur tiger. Front Microbiol. 2020;11:1665.
Oliveira BCM, Murray M, Tseng F, Widmer G. The fecal microbiota of wild and captive raptors. Anim Microbiome. 2020;2(1):15.
Prabhu VR, Wasimuddin, Kamalakkannan R, Arjun MS, Nagarajan M. Consequences of domestication on gut microbiome: a comparative study between wild gaur and domestic mithun. Front Microbiol. 2020;11:133.
Bowerman KL, Knowles SCL, Bradley JE, Baltrūnaitė L, Lynch MDJ, Jones KM, Hugenholtz P. Effects of laboratory domestication on the rodent gut microbiome. ISME Commun. 2021;1(1):1–14.
Liu C, Hu J, Wu Y, Irwin DM, Chen W, Zhang Z, Yu L. Comparative study of gut microbiota from captive and confiscated-rescued wild pangolins. J Genet Genom. 2021;48(9):825–35.
Sawaswong V, Praianantathavorn K, Chanchaem P, Khamwut A, Kemthong T, Hamada Y, Malaivijitnond S, Payungporn S. Comparative analysis of oral-gut microbiota between captive and wild long-tailed macaque in Thailand. Sci Rep. 2021;11(1):14280.
Carey HV, Walters WA, Knight R. Seasonal restructuring of the ground squirrel gut microbiota over the annual hibernation cycle. Am J Physiol Regul Integr Comp Physiol. 2013;304(1):R33-42.
Dill-McFarland KA, Neil KL, Zeng A, Sprenger RJ, Kurtz CC, Suen G, Carey HV. Hibernation alters the diversity and composition of mucosa-associated bacteria while enhancing antimicrobial defence in the gut of 13-lined ground squirrels. Mol Ecol. 2014;23(18):4658–69.
Stevenson TJ, Buck CL, Duddleston KN. Temporal dynamics of the cecal gut microbiota of juvenile arctic ground squirrels: a strong litter effect across the first active season. Appl Environ Microbiol. 2014;80(14):4260–8.
Stevenson TJ, Duddleston KN, Buck CL. Effects of season and host physiological state on the diversity, density, and activity of the arctic ground squirrel cecal microbiota. Appl Environ Microbiol. 2014;80(18):5611–22.
Sonoyama K, Fujiwara R, Takemura N, Ogasawara T, Watanabe J, Ito H, Morita T. Response of gut microbiota to fasting and hibernation in Syrian hamsters. Appl Environ Microbiol. 2009;75(20):6451–6.
Luo D, Ziebell S, An L. An informative approach on differential abundance analysis for time-course metagenomic sequencing data. Bioinformatics. 2017;33(9):1286–92.
Nelson MC, Morrison HG, Benjamino J, Grim SL, Graf J. Analysis, optimization and verification of Illumina-generated 16S rRNA gene amplicon surveys. PLoS ONE. 2014;9(4):e94249.
Sinclair L, Osman OA, Bertilsson S, Eiler A. Microbial community composition and diversity via 16S rRNA gene amplicons: evaluating the illumina platform. PLoS ONE. 2015;10(2):e0116955.
Good IJ. The population frequencies of species and the estimation of population parameters. Biometrika. 1953;40(3/4):237–64.
La Rosa PS, Brooks JP, Deych E, Boone EL, Edwards DJ, Wang Q, Sodergren E, Weinstock G, Shannon WD. Hypothesis testing and power calculations for taxonomic-based human microbiome data. PLoS ONE. 2012;7(12):e52078.
Xia Y, Sun J, Ding-Geng C. Power and sample size calculations for microbiome data. In: Chen J, Ding-Geng C, editors. Statistical analysis of microbiome data with R. Singapore: Springer; 2018.
Song H, Kim J, Guk JH, Kim WH, Nam H, Suh JG, Seong JK, Cho S. Metagenomic analysis of the gut microbiota of wild mice, a newly identified reservoir of Campylobacter. Front Cell Infect Microbiol. 2020;10:596149.
Derrien M, van Passel MWJ, van de Bovenkamp JHB, Schipper RG, de Vos WM, Dekker J. Mucin-bacterial interactions in the human oral cavity and digestive tract. Gut Microbes. 2010;1(4):254–68.
Ju T, Kong JY, Stothard P, Willing BP. Defining the role of Parasutterella, a previously uncharacterized member of the core gut microbiota. ISME J. 2019;13(6):1520–34.
Di Rienzi SC, Sharon I, Wrighton KC, Koren O, Hug LA, Thomas BC, Goodrich JK, Bell JT, Spector TD, Banfield JF, et al. The human gut and groundwater harbor non-photosynthetic bacteria belonging to a new candidate phylum sibling to Cyanobacteria. Elife. 2013;2:e01102.
Soo RM, Skennerton CT, Sekiguchi Y, Imelfort M, Paech SJ, Dennis PG, Steen JA, Parks DH, Tyson GW, Hugenholtz P. Photosynthesis is not a universal feature of the phylum Cyanobacteria. PeerJ. 2014. https://doi.org/10.7287/peerj.preprints.204v2.
Utami YD, Kuwahara H, Murakami T, Morikawa T, Sugaya K, Kihara K, Yuki M, Lo N, Deevong P, Hasin S, et al. Phylogenetic diversity and single-cell genome analysis of “Melainabacteria”, a non-photosynthetic cyanobacterial group, in the termite gut. Microbes Environ. 2018;33(1):50–7.
Flint HJ, Scott KP, Duncan SH, Louis P, Forano E. Microbial degradation of complex carbohydrates in the gut. Gut Microbes. 2012;3(4):289–306.
Kaoutari AE, Armougom F, Gordon JI, Raoult D, Henrissat B. The abundance and variety of carbohydrate-active enzymes in the human gut microbiota. Nat Rev Microbiol. 2013;11(7):497–504.
Pengelley ET, Fisher KC. Rhythmical arousal from hibernation in the golden-mantled ground squirrel, Citellus lateralis tescorum. Can J Zool. 1961;39:105–20.
Yuen KK, Dixon WJ. The approximate behaviour and performance of the two-sample trimmed t. Biometrika. 1973;60(2):369.
Yuen KK. The two-sample trimmed t for unequal population variances. Biometrika. 1974;61(1):165.
Skarlupka JH, Kamenetsky ME, Jewell KA, Suen G. The ruminal bacterial community in lactating dairy cows has limited variation on a day-to-day basis. J Anim Sci Biotechnol. 2019;10:66.
Streubel DP, Fitzgerald JP. Spermophilus tridecemlineatus. Mamm Species. 1978;103:1–5.
Stevenson DM, Weimer PJ. Dominance of Prevotella and low abundance of classical ruminal bacterial species in the bovine rumen revealed by relative quantification real-time PCR. Appl Microbiol Biotechnol. 2007;75(1):165–74.
Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Lozupone CA, Turnbaugh PJ, Fierer N, Knight R. Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. PNAS. 2011;108(1):4516–22.
Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD. Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Appl Environ Microbiol. 2013;79(17):5112–20.
Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, Lesniewski RA, Oakley BB, Parks DH, Robinson CJ, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75(23):7537–41.
Pruesse E, Quast C, Knittel K, Fuchs BM, Ludwig W, Peplies J, Glockner FO. SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids Res. 2007;35(21):7188–96.
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, Peplies J, Glockner FO. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41(Database issue):D590-596.
Yilmaz P, Parfrey LW, Yarza P, Gerken J, Pruesse E, Quast C, Schweer T, Peplies J, Ludwig W, Glockner FO. The SILVA and “All-species Living Tree Project (LTP)” taxonomic frameworks. Nucleic Acids Res. 2014;42(Database issue):D643-648.
Glöckner FO, Yilmaz P, Quast C, Gerken J, Beccati A, Ciuprina A, Bruns G, Yarza P, Peplies J, Westram R, et al. 25 Years of serving the community with ribosomal RNA gene reference databases and tools. J Biotechnol. 2017;261:169–76.
R Development Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2020.
McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE. 2013;8(4):e61217.
Revell LJ. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol Evol. 2012;3(2):217–23.
Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer-Verlag; 2009.
Auguie B. gridExtra: miscellaneous functions for "Grid" graphics, 2.3 ed. 2017.
Paradis E, Schliep K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics. 2019;35(3):526–8.
Wickham H, François R, Henry L, Müller K. dplyr: a grammar of data manipulation, 0.8.0.1 edn. 2019.
Dinno A: dunn.test: Dunn's test of multiple comparisons using rank sums, 1.3.5 ed. 2017.
Torchiano M. effsize: efficient effect size computation. R package version 0.8.1 ed. 2020.
Robinson A. equivalence: provides tests and graphics for assessing tests of equivalence, 0.7.2 ed. 2016.
Kembel SW, Cowan PD, Helmus MR, Cornwell WK, Morlon H, Ackerly DD, Blomberg SP, Webb CO. Picante: R tools for integrating phylogenies and ecology. Bioinformatics. 2010;26(11):1463–4.
Champely S. pwr: basic functions for power analysis. R package version 1.3–0 edn. 2020.
Wickham H, tidyr: tidy messy data, 1.1.2 ed. 2020.
Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O'Hara RB, Simpson GL, Solymos P, et al. vegan: community ecology package, 2.4–6 ed. 2018.
Shannon CE. A mathematical theory of communication. Bell Syst Tech J. 1948;27(379–423):623–56.
Tucker CM, Cadotte MW, Carvalho SB, Davies TJ, Ferrier S, Fritz SA, Grenyer R, Helmus MR, Jin LS, Mooers AO, et al. A guide to phylogenetic metrics for conservation, community ecology and macroecology. Biol Rev Camb Philos Soc. 2017;92(2):698–715.
Anderson MJ, Walsh DCI. PERMANOVA, ANOSIM, and the Mantel test in the face of heterogenous dispersions: what null hypothesis are you testing? Ecol Monogr. 2013;83(4):557–74.
The authors thank Michael Grahn for technical assistance with animals, Dr. Kimberly A. Dill-McFarland and Andrew J. Steinberger for technical assistance with sample collection, and Madison S. Cox and Joseph H. Skarlupka for their PCR troubleshooting wisdom. The authors also thank the Suen Lab for their support and comments on the manuscript.
This work was supported by National Science Foundation Grant IOS-1558044 to HVC and GS. EC was supported by the National Institute of General Medical Sciences of the National Institutes of Health traineeship under award number T32GM008349, the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1747503, and the Graduate School and the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin-Madison with funding from the Wisconsin Alumni Research Foundation.
All procedures were approved by the University of Wisconsin-Madison Institutional Animal Care and Use Committee under protocols V001229 and V005481.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Alpha diversity of content and mucosa microbiotas within each experimental group. Violin plots display the number of observed OTUs (first column), Faith’s phylogenetic diversity (second column), phylogenetic evenness (mean pairwise distance (MPD); third column), and Shannon’s diversity (fourth column) for each experimental group: (A – D) Summer Wild, (E – H) Summer Captive, (I – L) Torpor, (M – P) IBA, (Q – T) Spring. There were no significant comparisons (all t-test adj P ≥ 0.889). Figure S2. Beta diversity of content and mucosa microbiotas within each experimental group. Weighted UniFrac (first column), unweighted UniFrac (second column), and Bray-Curtis dissimilarity (third column) are displayed on PCoA ordinations for each group: (A – C) Summer Wild, (D – F) Summer Captive, (G – I) Torpor, (J – L) IBA, (M – O) Spring. Figure S3. Alpha diversity between summer Captive and Wild content microbiotas. Violin plots display four alpha diversity metrics: (A) the number of observed OTUs, (B) Faith’s phylogenetic diversity, (C) phylogenetic evenness (MPD), and (D) Shannon’s diversity. An asterisk indicates a significant difference (t-test adj P < 0.05) and “ns” indicates no significant difference (adj P > 0.05). Figure S4. Beta diversity of summer Captive and Wild content microbiotas. (A) Weighted UniFrac, (B) unweighted UniFrac, and (C) Bray-Curtis dissimilarity are displayed on principal coordinate analysis PCoA ordinations. Groups are depicted with different colors. Figure S5. Alpha diversity comparison among content microbiotas across captive groups. Violin plots display four alpha diversity metrics: (A) the number of observed OTUs, (B) Faith’s phylogenetic diversity, (C) phylogenetic evenness (MPD), and (D) Shannon’s diversity. Groups that share a letter are not significantly different (Tukey’s HSD adj P > 0.05), whereas groups that share no letters are significant different (adj P < 0.05). Metrics with no significant comparisons between groups are indicated with “ns." Figure S6. Beta diversity of content microbiotas in captive groups. (A) Weighted UniFrac, (B) unweighted UniFrac, and (C) Bray-Curtis dissimilarity are displayed on principal coordinate analysis PCoA ordinations. Groups are depicted with different colors.
Table S1. Squirrel metadata. Table S2. Gut microbiota comparison of squirrels that received isoflurane or CO2 for anesthesia. To determine whether anesthesia had a significant impact on the gut microbiota, we examined IBA squirrels as this was the only group in which more than one anesthesia method was used. We compared IBA animals that received isoflurane or CO2 using the robust two one-sided test (RTOST) of equivalence [84, 102,103,104]. RTOST is a nonparametric test whose null and alternative hypotheses are switched compared to that of a traditional t-test: the null hypothesis is that two means are not equivalent and the alternative hypothesis is that two means are equivalent. Using OTU tables, each sample from a squirrel that received CO2 was compared to the mean of all samples from squirrels that received isoflurane (Table S1). P-values were adjusted for false discovery rate using the Benjamini-Hochberg procedure. All comparisons were significant; therefore, we rejected our null hypothesis that the means (or microbiotas of isoflurane- and CO2 - administered squirrels) are not equivalent. Table S3. Metadata and alpha diversity metrics for the four outlier samples. Table S4. Effect size and power for alpha diversity comparisons between paired content and mucosa samples in each experimental group. Effect sizes were calculated in R using cohen.d (effsize package ), and power was calculated using pwr.t2n.test (pwr package ). Table S5. Effect size and power for Wild vs. Captive group comparisons of alpha diversity and phyla relative abundances. Effect sizes were calculated in R using cohen.d (effsize package ), and power was calculated using pwr.t2n.test (pwr package ). Table S6. PERMANOVA results from beta diversity comparison of content and mucosa beta diversity in each experimental group. Three beta diversity metrics (weighted UniFrac, unweighted UniFrac, and Bray-Curtis dissimilarity) and two variables (squirrel ID and sample type) were tested. Table S7. Phyla relative abundances in summer Captive and Wild content microbiotas. Phyla have total relative abundances > 1% and are displayed from overall highest relative abundance to lowest relative abundance. Relative abundances are displayed as mean ± standard error and adjusted p-values are from t-tests or Wilcoxon Rank Sum tests, depending on data normality. Table S8. OTUs that significantly contribute to Bray-Curtis dissimilarity between Captive and Wild content microbiotas. OTUs were identified with SIMPER and statistically tested with Kruskal-Wallis tests. P-values were corrected for false discovery rate using the Benjamini-Hochberg procedure. Only OTUs that accounted for ≥ 1% of the differences in beta diversity and had adj P < 0.05 were considered significant. Table S9. Core OTUs in summer Captive and Wild mucosa microbiotas. X's in a Captive or Wild column indicate that the OTU is a core OTU in that group, while empty cells indicate that the OTU is not a core OTU in the group. Table S10. Core OTUs in summer Captive and Wild content microbiotas. X's in a Captive or Wild column indicate that the OTU is a core OTU in that group, while empty cells indicate that the OTU is not a core OTU in the group. Table S11. Phyla relative abundances in content microbiotas of captive groups. Phyla have total relative abundance >1% and are displayed from overall highest relative abundance to lowest relative abundance. Relative abundances are displayed as mean ± standard error and adjusted p-values are from ANOVA and Tukey's HSD tests, or Kruskal-Wallis and Dunn's tests, depending on data normality. Table S12. Core OTUs in mucosa microbiotas across captive groups. X's in a captive group column indicate that the OTU is a core OTU in that group, while empty cells indicate that the OTU is not a core OTU in the group. Table S13. Core OTUs in content microbiotas across captive groups. X's in a captive group column indicate that the OTU is a core OTU in that group, while empty cells indicate that the OTU is not a core OTU in the group.
About this article
Cite this article
Chiang, E., Deblois, C.L., Carey, H.V. et al. Characterization of captive and wild 13-lined ground squirrel cecal microbiotas using Illumina-based sequencing. anim microbiome 4, 1 (2022). https://doi.org/10.1186/s42523-021-00154-9