Animals
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. 1). 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:12hr 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.
Sampling Regimen
The 26 captive TLGS were randomly assigned to one of four groups (Fig. 1): 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 [31]. 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 Table S1.
Sample collection
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 (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 (Table S2). We did not find evidence that the microbiotas of isoflurane- or CO2-administered animals were different. The methods and results for this analysis are shown in 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 [32]. 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.
DNA Extraction
Total genomic DNA was extracted from each sample using a phenol:chloroform extraction protocol [33] 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 (forward: GTGCCAGCMGCCGCGGTAA, reverse: GGACTACHVGGGTWTCTAAT) [34]. 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 x 250 v2 kit.
Microbiota sequence clean-up
Sequences were demultiplexed on the Illumina MiSeq and processed using mothur v1.44.3 [35] and SILVA release 132 [36–39]. 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/mothur_outputs.
We followed the mothur SOP accessed on March 19, 2021 with the following modifications: we used screen.seqs maxlength = 300 and chimera.uchime. Sequences were grouped into 97% operational taxonomic units (OTUs). Samples were normalized to 11,013 sequences. This cutoff was selected based on the sample with the lowest number of sequences that resulted in all normalized samples having Good’s coverage ≥ 95%, which indicated sufficient sequencing [40]. 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.
Preparing data for statistical analysis in R
Outputs from mothur were imported into R version 3.4.3 [41] using the phyloseq [42] and phytools [43] packages. All visualizations were created using the ggplot2 [44], grid [41], and gridExtra [45] packages. The following packages were also used for statistical analyses: ape [46], dplyr [47], dunn.test [48], equivalence [49], picante [50], stats [41], tidyr [51], and vegan [52].
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. 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 unique OTUs, Faith’s phylogenetic diversity, phylogenetic evenness, and Shannon’s weighted diversity index [53]. The number of unique OTUs was calculated using the phyloseq package [42]. Faith’s phylogenetic diversity was computed by taking the sum of all branches in a tree using the picante package (pd) [50]. Because phylogenetic diversity is positively correlated with richness, we included a metric unbiased by richness: mean pairwise distance (MPD), which measures phylogenetic evenness [54]. 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 [41]; ses.mpd, picante package [50]). Lastly, Shannon’s weighted diversity index was calculated in mothur [35]. We tested the normality of each alpha diversity metric using quantile-quantile plots and the Shapiro test (qqnorm, qqline, shapiro.test, stats package [41]) 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 [41]) followed by pairwise Tukey’s Honest Significant Difference tests (TukeyHSD, stats package [41]). Mucosa MPD and Shannon’s diversity index were also analyzed this way. Mucosa number of unique OTUs and phylogenetic diversity were not normally distributed and were therefore analyzed using Kruskal-Wallis (kruskal.test, stats package [41]) followed by Dunn’s Test (dunn.test, dunn.test package [48]). All p-values were adjusted for false discovery rate using the Benjamini-Hochberg procedure (p.adjust, stats package [41]). 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) [42] and visualized using principal coordinate analysis (PCoA) ordinations (ordinate, phyloseq package [42]). To test whether there were significant differences in variance/dispersion, homogeneity of groups dispersions tests (betadisper, vegan package [52]) were used. Next, we tested if beta diversity centroids differed between experimental groups using permutational multivariate analysis of variance [55] (PERMANOVA; adonis, vegan package [52]). 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 [52]) 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 analysis
Phylum-level relative abundances were calculated using the phyloseq [42] and dplyr [47] 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 [41]). 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 (Table S4).
Identify Core OTUs
Core OTUs were defined as OTUs that were present in every sample within a group and were identified using the phyloseq package [42] (Table S8). Core OTUs shared among groups were identified by taking the intersect of the core OTUs in two or more groups.