Microbial community structure and composition is associated with host species and sex in Sigmodon cotton rats

Background The cotton rat (genus Sigmodon) is an essential small animal model for the study of human infectious disease and viral therapeutic development. However, the impact of the host microbiome on infection outcomes has not been explored in this model, partly due to the lack of a comprehensive characterization of microbial communities across different cotton rat species. Understanding the dynamics of their microbiome could significantly help to better understand its role when modeling viral infections in this animal model. Results We examined the bacterial communities of the gut and three external sites (skin, ear, and nose) of two inbred species of cotton rats commonly used in research (S. hispidus and S. fulviventer) by using 16S rRNA gene sequencing, constituting the first comprehensive characterization of the cotton rat microbiome. We showed that S. fulviventer maintained higher alpha diversity and richness than S. hispidus at external sites (skin, ear, nose), but there were no differentially abundant genera. However, S. fulviventer and S. hispidus had distinct fecal microbiomes composed of several significantly differentially abundant genera. Whole metagenomic shotgun sequencing of fecal samples identified species-level differences between S. hispidus and S. fulviventer, as well as different metabolic pathway functions as a result of differential host microbiome contributions. Furthermore, the microbiome composition of the external sites showed significant sex-based differences while fecal communities were not largely different. Conclusions Our study shows that host genetic background potentially exerts homeostatic pressures, resulting in distinct microbiomes for two different inbred cotton rat species. Because of the numerous studies that have uncovered strong relationships between host microbiome, viral infection outcomes, and immune responses, our findings represent a strong contribution for understanding the impact of different microbial communities on viral pathogenesis. Furthermore, we provide novel cotton rat microbiome data as a springboard to uncover the full therapeutic potential of the microbiome against viral infections. Supplementary Information The online version contains supplementary material available at 10.1186/s42523-021-00090-8.


Background
The commensal microbiome can dramatically in uence many aspects of host health and disease, such as homeostatic signaling, nutrient acquisition, and protection from or exacerbation of infections (1)(2)(3).
The majority of early studies established that environmental factors play the major role in shaping and modulating the host microbiome (4)(5)(6)(7). In addition, recent studies have emphasized that there is a signi cant role of host genetics on co-evolution of the host and its associated microbiome (8)(9)(10). For example, the murine genetic background is a stronger determinant of microbiome composition and structure than environmental stimuli (11). Similarly, genetic polymorphisms, heritability, and overall host genetics in humans can also shape how commensal bacteria evolve alongside the host (12)(13)(14). The microbiome has also been instrumental in predicting and protecting against severe viral disease outcomes (15,16). However, this burgeoning eld of bacteria-host-virus interactions has been limited by a lack of translational models to study mechanisms of virus-microbiome interaction.
Cotton rats (genus Sigmodon) are an important small animal model to study various respiratory diseases, including Respiratory Syncytial Virus (RSV) (17), in uenza A Virus (IAV) (18,19), parain uenza virus (20,21), measles (22), human metapneumovirus (23), enterovirus (24), and rhinovirus (25) due to comparable human disease outcomes (26). Cotton rats have also provided a useful model for nasal colonization studies (especially with Staphylococcus aureus) due to their human-like nasal histology (27). Furthermore, cotton rats are a useful tool for research since they harbor zoonotic viruses in the wild like Alphavirus (Equine Encephalitis virus), Hantavirus (Black Creek Canal virus, Bayou virus), Cardiovirus, Arenavirus (Tamiami virus), and Flavivirus (West Nile virus) (28)(29)(30)(31)(32)(33). However, viral-microbiota interactions in this model have not been explored due to the lack of a comprehensive analysis of the normal microbiota in healthy cotton rats and microbiome comparison between different species used in laboratory settings. To date, only one study has examined the nasal micro ora of healthy S. hispidus but was limited by the sample number and lack of longitudinal timepoints (34).
To comprehensively characterize and establish the structure and composition of the cotton rat microbiome, we collected longitudinal samples from four different body sites of two commonly used inbred cotton rat species, Sigmodon hispidus and S. fulviventer, maintained under the same environment and dietary conditions. Our microbiome characterization using both 16S rRNA gene and whole metagenomic sequencing (WMS) comprehensively establishes the microbiome community structure and composition of different body sites in cotton rats and showed distinct structure based on the cotton rat species and sex. WMS also showed differential metabolic potential of the community between species.
Overall, this study not only adds to the small but rapidly expanding literature of the in uence of host genetics on the microbiome, but also describes an appropriate animal model for studying microbiome in uences on viral and bacterial diseases.

Results
Characterization of cotton rat microbiome from multiple body sites.
Two groups of 10 young male cotton rats of S. fulviventer and S. hispidus were observed longitudinally for 111 days to characterize the healthy cotton rat microbiome structure and composition. A total of 140 samples were collected and processed for 16S rRNA gene sequencing: ear swabs (20 swabs/day 95), nasal brushes (20 swabs/day 95), skin swabs (20 swabs/day 95), and fecal samples (80 swabs/days 0, 4, 34, and 111) (Fig. 1A). DNA was extracted, and the V4 region of the 16S rRNA gene was ampli ed then sequenced on the Illumina MiSeq platform with 2 × 250 base pair reads, generating an average of 35,194 reads per sample with 96.4% of samples passing quality control (i.e., > 1000 reads/sample) for analysis.
Microbiome data was processed by following the mothur SOP, and operational taxonomic units (OTUs) were clustered at 97% identity. We then examined the association of the alpha diversity, beta diversity, and abundance of taxa at each site using R version 3.5.1.
Differences in the microbiome community structure and composition between cotton rat species. Phyla within S. fulviventer and S. hispidus external sites (i.e., ear, nose, skin) were similarly dominated by Proteobacteria, Actinobacteria, Tenericutes, and Firmicutes. Fecal communities consisted mostly of two dominant phyla with opposite abundances between cotton rat species: higher Bacteroidetes abundance in S. fulviventer compared to S. hispidus (50.0% vs. 42.4% respectively) and higher Firmicutes abundance in S. hispidus compared to S. fulviventer (36.2% vs. 31.8% respectively) ( Table 1). The distribution of top 20 genera in each cotton rat species is shown in Figure S1.
We found that the cotton rat gut microbiome was stable over time in both cotton rat species. Richness and alpha diversity did not signi cantly change over time, no taxa were signi cantly differentially abundant when experimental day was set as the outcome variable, and beta diversity testing revealed no signi cant shifts in microbiome composition over time (data not shown). Subsequently, we analyzed groups by computing mean counts for individual cotton rats across time points. Comparison of individual body site microbiomes of both S. fulviventer and S. hispidus across all time points showed community distinctions between species. All sites from S. fulviventer consistently had higher richness (Observed OTUs, S.chao1 index) and alpha diversity (Shannon Index, Simpson's Index) when compared to S. hispidus ( Fig. 1B and C, Figure S2; all values p < 0.05). Beta diversity, computed by calculating Bray-Curtis dissimilarity between samples at the OTU level, showed unique composition between the fecal communities of S. fulviventer and S. hispidus ( Fig. 1D; PERMANOVA p = 0.00025, permutations = 4000, R 2 = 0.26284; beta dispersion, p = 0.00025, R 2 = 0.34564, Figure S3, S4D). However, comparison of beta diversity metrics of individual external sites from S. fulviventer and S. hispidus did not show signi cant differences ( Fig. 1D; Figure S4A-C).
Escherichia/Shigella was more abundant in the S. hispidus gut (log 2 fold change = 7.30, q = 3.79E-11) but had very low relative abundance compared to other genera. S. fulviventer had a higher abundance of 18 unique genera in the gut (q < 0.05), including Clostridium sensu stricto (log 2 fold change=-7.19, q = 7.77E-12), Elusimicrobium (log 2 fold change=-6.22, q = 8.04E-18), and Hespellia (log 2 fold change=-5.11, q = 2.21E-04) ( Fig. 2A). Full data of differentially abundant taxa at both the genus and family levels are shown in Supplemental File 1. A total of 32 of the 36 DESeq2-calculated differentially abundant genera were also con rmed using the GeneSelector (36) R package ( Figure S5, Supplemental File 2). A stability selection model showed Lactobacillus as one of the top genera (as well as those unclassi ed within the phyla Bacteroidetes) with a high probability of predicting whether a fecal sample was from S. hispidus or S. fulviventer (Fig. 2B). While Lactobacillus was one of the top 20 most abundant bacteria of the skin, ear, and nose microbiomes of both S. fulviventer and S. hispidus ( Figure S1), it was not signi cantly differentially abundant between the two species at any other sites except feces (Fig. 2C).
There were few speci c taxa at both the genus and family levels that were signi cantly differentially abundant (p < 0.05, q < 0.10) between cotton rat species when examining external communities (ear, nose, skin). Enterobacteriaceae was more abundant in the S. hispidus nose (log 2 fold change = 1.21, q = 1.05E-06), ear (log 2 fold change = 0.48, q = 0.022), and skin (log 2 fold change = 0.47, q = 0.0031) compared to the equivalent sites in S. fulviventer. Corynebacteriaceae was more abundant in the S. hispidus nose (log 2 fold change = 0.532, q = 0.0305) and ear (log 2 fold change = 0.571, q = 0.0470) compared to the equivalent sites in S. fulviventer. Leptotrichiaceae was more abundant in the S. fulviventer nose (log 2 fold change=-1.41, q = 0.0456) and ear (log 2 fold change=-1.41, q = 0.047) compared to the equivalent sites in S. hispidus. Pseudomonas was more abundant in the S. hispidus ear (log 2 fold change = 1.33, q = 0.070) than in S. fulviventer. Barnesiella and Porphyromonadaceae were more abundant in the S. fulviventer ear (log 2 fold change=-4.58, q = 0.066; log 2 fold change=-2.11, q = 0.063). Leuconostocaceae was more abundant in the S. fulviventer nose (log 2 fold change=-1.56, q = 0.064). Full data at both genus and family levels are shown in Supplementary File 1.
Con rmation of 16S rRNA gene sequencing data using traditional culture methods.
For quantitative comparison of bacterial load between cotton rat species, we used qPCR analysis of total bacterial DNA extracted from homogenized stool (equal weight/volume) and found that the bacterial load was signi cantly higher in S. hispidus than S. fulviventer (Fig. 2C). We also plated an aliquot of normalized, homogenized stool on Lactobacillus-speci c agar (De Man, Rogosa and Sharpe agar) and observed that the number of colony-forming units (CFU) per gram in S. hispidus stool was signi cantly higher than in S. fulviventer stool (Fig. 2D). We found that 86% of colonies grown from S. hispidus stool were Lactobacillus-positive, compared to zero Lactobacillus-positive colonies grown from S. fulviventer stool (Fig. 2E). Sequencing of colonies showed Lactobacillus gasseri and Lactobacillus. reuteri to be the two prominent bacterial species found in S. hispidus stool (Fig. 2E). This signi cant trend supports the relative abundance of Lactobacillus as determined by 16S rRNA gene sequencing, where Lactobacillus was signi cantly more abundant in S. hispidus compared to S. fulviventer (Fig. 2F). Additionally, Corynebacterium and Bacteroides species were also identi ed in S. hispidus stool samples. Sanger sequencing of isolates from S. fulviventer stool identi ed the presence of Enterococcus gallinarum and E. casselifavus.
Differences in the microbiome community structure and composition based on sex.
We conducted a follow-up experiment to assess if there were any associations between host sex and microbiome community structure and composition. This cohort included 13 S. fulviventer cotton rats (10 males, 3 females) and 16 S. hispidus cotton rats (5 males, 9 females). Animals in both groups were 4-6 weeks old and weighed approximately 100 g and were observed for 28 days, with fecal samples collected on days 0, 7, 13, 21, and 28 and nose, ear, and skin swabs collected on days 7 and 28. We performed 16S rRNA gene sequencing to examine any effect of host sex. Alpha diversity metrics indicated signi cant differences in richness (Observed OTUs, Chao1) and diversity (Shannon Index, Simpson Index) between male and female S. fulviventer at both the ear and fecal microbiomes ( Figure S6A-D), but there were no signi cant differences in richness and diversity of the microbiomes of male and female cotton rats in the nose and skin for both S. fulviventer and S. hispidus (with the exception of S. hispidus skin diversity, Figure S6D). Overall, differences between host sex were most pronounced in the gut compared to external sites ( Figure S6) but only in S. fulviventer. Beta-diversity measurements of each species revealed that microbial composition of the gut was signi cantly dissimilar between male and female cotton rats for both S. fulviventer and S. hispidus (Fig. 3A, B; S. hispidus PERMANOVA p = 0.00025, beta-dispersion p = 0.1116; S. fulviventer PERMANOVA p = 0.00025, beta-dispersion p = 7E-04). There were also notable differences between sex at S. hispidus skin and nose (Fig. 3E, G). Differential abundance analysis using DESeq2 was conducted between males and females at each site (Table S1). While the fecal community structure differed, there were only 3 differentially abundant genera due to sex in the S. hispidus gut and no different genera in S. fulviventer. There were differential taxa between sexes at external sites of both S. hispidus (21 genera) and S. fulviventer (13 genera). Full results at genus and family levels are listed in Supplementary File 3.
Differences in the microbiome between cotton rat species assessed by whole metagenomic sequencing.
Due to the dramatic differences between the S. fulviventer and S. hispidus gut microbiomes detected by 16S rRNA gene sequencing, we pursued further analysis in order to understand community differences at the species and strain level as well as differences in microbiome functional potential. DNA extracted from 10 male cotton rat stool samples (S. hispidus = 5, S. fulviventer = 5) at both days 34 and 111 (20 samples total) from the rst experimental group were processed for shotgun metagenomic sequencing, which generated 1.11 × 10 9 reads (5.57 × 10 7 average reads per sample), comprising of 334,037 megabases (16,701 average megabases per sample) and 17.2% duplicate reads.
Whole metagenomic sequencing data showed differences at the species level that validated the 16S rRNA gene sequencing data. Abundances of several bacterial species were found to be statistically different (q < 0.05) between cotton rat species based on taxonomic classi cation as performed by MetaPhlAn2 followed by differential abundance analysis by both hierarchical clustering (based on Bray-Curtis dissimilarity) of the top 25 most abundant species (Fig. 4A) and DESeq2 (Supplemental File 4). Lactobacillus reuteri, L. gasseri, and the novel L. sp. ASF360 predominated the gut of S. hipsidus (Fig. 4A), and many other Lactobacillus species were signi cantly more abundant in S. hispidus samples compared to S. fulviventer ( Figure S7). Total Lactobacillus within the S. fulviventer gut was signi cantly less abundant but included species unique to S. fulviventer, including L. murinus, L. BHWM-4, and L. animalis. Akkermansia muciniphilia was signi cantly more abundant in S. fulviventer compared to S. hispidus. Ruminococcus torques, Helicobacter cinaedi, and Oscillibacter spp. were of higher abundance in the S. hispidus gut. Parabacteroides spp. (including P. johnsonii) and Odoribacter were more abundant in the S. fulviventer gut (Supplementary File 4). Proportional counts and raw counts can be found in Supplementary Files 5 and 6 respectively. Differential functional potential between cotton rat species microbiome.
To understand the biological implications of these differences, Humann2 (37) was used to map any functional differences (MetaCyc pathway database (38)) de ned by identi ed gene families and bacterial pro les. We identi ed 415 pathways (with nearly all associated with bacteria present in the sample) in the two cotton rat species based on the MetaCyc database (Fig. 4B). The majority of these pathways included biosynthesis (39.60%) and degradation/utilization/assimilation (18.79%) pathways, as well as several overarching superpathways (27.18%) and energy/metabolite production pathways (10.57%). More speci cally, these pathways were part of several instrumental superclass ontologies that metabolize (including de novo pathways) electron carriers, vitamins, fatty acids, lipids, amino acids, carbohydrates, secondary metabolites, and fermentation-derived energy (Fig. 4C). Interestingly, several pathways were differentially abundant between cotton rat species. Each cotton rat species had unique pathways contributed to by their microbiomes (S. fulviventer = 14, and S. hispidus = 27), and most of these involved biosynthesis (Supplemental File 7).
In relation to differentially abundant bacteria species, we found that 44 pathways were solely driven by Lactobacillus gasseri, L. reuteri, and L. ASF360 by matching reads from MetaPhlAn2 bacterial identi cations with Humann2 predicted pathways. Several of these pathways were more highly expressed in S. hispidus (Fig. 5). These included L-proline biosynthesis from arginine (catalyzed by bacterial enzymes), inosine-5'-phosphate biosynthesis (for de novo synthesis of purines), pyruvate fermentation to acetate/lactate (for anaerobic energy production), adenosine deoxyribonuclease de novo biosynthesis (to promote ADP production), and D-galactose degradation (breakdown of D-galactose to a useable form in glycolysis). Akkermansia municiphilia was the driver of 25 other pathways, many of which were highly expressed in S. fulviventer compared to S. hispidus ( Figure S8). These included Lisoleucine biosynthesis (for production of leucine and isoleucine), phosphopantothenate biosynthesis (to produce vitamin B5 de novo, of which animals cannot produce, and to feed production of coenzyme A and acyl carrier protein), glycolysis (particularly the degradation of starches for reductants and energy for anabolic pathways), and L-valine biosynthesis. Statistical comparison of all pathways can be found in Supplemental File 8.

Discussion
Here we have comprehensively characterized the cotton rat microbiome and compared bacterial communities of two different species (S. hispidus and S. fulviventer) that were housed in the same facility with identical diets. From these analyses, we were able to uncover species-speci c differences of gut bacteria, even though the two species had the same diet over many generations and were housed in separate cages in the same room. Interestingly, their external microbiomes (ear, nose, skin) were remarkably similar based on beta-diversity testing and testing for differential bacterial taxa but signi cantly different based on alpha diversity and richness measurements (which mimicked that of fecal communities). This data further supports that, while environmental factors play a vital role in shaping microbiome structure and composition, underlying host genetics exerts homeostatic pressure for distinct microbiomes between populations.
The most recent phylogenetic analysis of the genus Sigmodon sp. found that S. hispidus and S. fulviventer diverged 5.4 million years ago (39). In the wild, S. hispidus and S. fulviventer are sympatric species, with S. fulviventer being the more dominant animal (40). Separate inbreeding of the two species has made them a useful small animal model in laboratory research (26,41). Our data show that even when inbred and adapted to a controlled laboratory environment, each species still maintains a unique gut microbiome community structure and composition. This is important to note because the cotton rat is a useful model for respiratory viral infection, pathogenesis, and immunity research (26). Additionally, cotton rats are an excellent small animal model for pre-clinical studies for antiviral and vaccine candidates for viruses e.g., RSV (17), IAV (18,19), and human rhinovirus (HRV) (26). Recent studies have shown that gut microbiome composition and disruptions by antibiotics may affect some immune responses to u vaccination (42,43). Similarly, there are important correlations between vaccination status and the microbiome of the gut and lung (44,45). Further, it has also been shown that transplantation of the gut microbiome of wild-caught mice to laboratory mice can promote host tness and improve disease resistance and pathogenesis (46). In light of the above literature suggesting that the microbiome may play a key role in respiratory viral disease exacerbation or remediation and vaccine response (42)(43)(44)(45)47), our ndings exemplify the usefulness of this small animal model for understanding viral pathogenesis in the context of different microbial communities. Furthermore, the cotton rat model could be an optimal subject to uncovering and examining the full therapeutic potential of the microbiome by understanding how the host regulates and modulates bacterial communities.
While other small animal models including mice and rats have been used to explore the relationship between genetic patterns and bacterial homeostasis (11,(48)(49)(50), the cotton rat could be used in studying the interplay between differing host microbiomes and host immune responses to viral infections. For example, S. hispidus had a signi cantly higher amount of probiotic gut bacteria genera (Lactobacillus, Bi dobacterium) that have been associated with protection against the severe outcomes from RSV, IAV, and HRV (51)(52)(53). Other differentially abundant bacteria between the two species, such as Escherichia coli and Bacillus cereus, have also been shown to enhance both poliovirus and reovirus replication and pathogenesis (54). With the presence/absence of these bacteria in this animal model, along with our recently published annotated transcriptome (55), host responses in light of the microbiome could be elucidated. The cotton rat could also be an optimal model for supplementation studies of particular taxa in relation to these viruses.
Here, we have not only comprehensively characterized and established the key differences in microbiome community structure and function between multiple sites from two species of cotton rats housed in same facility for generations and fed with same diet, we also uncovered sex as a factor that can impact the microbiome composition of the gut. Our study has signi cant strengths compared to the lone study published so far (34), including a large sample size, longitudinal sampling, and shotgun metagenomic sequencing to characterize the gut microbiome at the species level. However, we acknowledge several shortcomings: 1) Due to the lack of an assembled cotton rat genome, we could not examine host genes, genetic patterns, or polymorphisms that may be driving differences in microbial colonization. 2) In both cohorts of cotton rats, samples were taken at different time points, and the cotton rats were followed for different durations. However, we found no evidence of changes in the gut microbiome across 111 days of sampling our rst cohort of adult cotton rats. 3) We only conducted shotgun metagenomic sequencing on the gut microbiome from male cotton rats from our rst cohort. While there is no cotton rat genome in the public databases, with the de novo assembly of the cotton rat transcriptome (55), further studies integrating the microbiome data and gene expression patterns may uncover more relevant information in regards to differences in host and microbiome interactions. We would also like to note that characterizing the microbiome of a unique animal species can be challenging due to the lack of host and/or microbiome sequence databases. The majority of databases and tools are that are commonly available are speci cally designed for human microbiome analyses and can result in misclassi cation of bacteria when used for new animal species. However, until we have better databases, interpretation of data has to be done with caution. Further research is warranted to understand species level microbiome differences and their impact on immune response in all small animal models for better interpretation of preclinical studies of vaccines and anti-microbiological agents. In spite of some limitations, our study creates a steppingstone for future research into these pressing questions of host-microbiome interactions during infection.

Conclusion
Overall, we have comprehensively characterized the cotton rat microbiome, an invaluable small animal model for viral and bacterial infections, and established key differences in microbiome community structure and function between multiple sites from two species of cotton rats (S. hispidus and S. fulviventer) housed in same facility for generations and fed with same diet. We also uncovered sex as a variable that can impact the microbiome composition of the gut. This foundational study establishes a foundation for future hypothesis testing experiments in understanding the role of microbiome in viral pathogenesis, especially for RSV and In uenza virus. Additionally, this study adds to the small but expanding literature in understanding the role of host genetics on microbiome composition and structure.

Animals
Four-to six-week-old cotton rats (~ 100 g) were obtained from the inbred colony maintained at Sigmovir Biosystems, Inc. (SBI). Cotton rats in the colony were seronegative by ELISA to adventitious respiratory viruses (i.e. Pneumonia Virus of Mice, Rat parvovirus, Rat coronavirus, Sendai virus). Animals were individually housed in large polycarbonate cages and fed a diet of standard rodent chow and water ad libitum.
For rigor and reproducibility, two independent animal experiments were carried out to characterize and establish the healthy cotton rat microbiome structure and composition by comparing two different species, S. fulviventer and S. hispidus. In the rst experimental group, 20 young male cotton rats were examined: S. fulviventer (n = 10) and S. hispidus (n = 10). Each animal was observed for 111 days, with nose, ear, and skin swabs collected at day 95 and fecal samples collected at days 0, 4, 34, and 111. These samples were used for microbiome characterization. To analyze any sex bias to the microbiome, a second experimental group (at a later time) included 13 young S. fulviventer (10 males, 3 females) and 16 young S. hispidus (5 males, 9 females). Healthy animals were monitored for 28 days with fecal samples collected on days 0, 7, 13, 21, and 28 and nose, ear, and skin swabs collected on days 7 and 28. To avoid ghting, all the animals were housed individually in large polycarbonate cages (with proper enrichment; nylon bone and glass jar). The cotton rat colony was maintained free of human and rodent viruses. All animal procedures followed NIH and USDA guidelines and were approved by the Sigmovir Biosystems, Inc. IACUC.

Sample Collection
Feces collection: One day prior to feces collection, cage beddings were changed in the late afternoon for each animal. Samples were collected between 10 am and 1 pm with sterile forceps. On average, 10-15 feces pellets were collected from each animal. Immediately after collection, samples were frozen at − 80 °C.
Nose swab: Sterile saline (~ 100 µl) was pipetted into both nostrils of anesthetized cotton rats positioned face down; Fisherbrand Sterile Swabs (Calcium Alginate Fiber Tipped Ultra ne Aluminum Applicator Swab) were then immediately placed in the nostrils to absorb the saline. Swabs were broken into sterile DNase/RNase-free 1.5 ml tubes and stored at − 80 °C.
Ear swab: Sterile saline (~ 100 µl) was pipetted up and down into both ears of each anesthetized cotton rat while the animal was kept in an anesthesia chamber for 1-2 minutes, and residual liquid was absorbed from each ear with Beaver Visitec Ultracell PVA Eye Spears pack of 5 (intended for uid absorption and tissue manipulation). Swabs were broken into sterile DNase/RNase-free 1.5 ml tubes and stored at − 80 °C.
Skin swab: Sterile saline (~ 200 µl) was put at the back of each anesthetized cotton rat (at approximately the same site for each animal) and rubbed vigorously using Fisherbrand Sterile Swabs (Calcium Alginate Fiber Tipped wood applicator swab). Swabs were broken in a sterile DNase/RNase-free 1.5 ml tubes and stored at − 80 °C.
Microbiome DNA extraction and 16S rRNA gene sequencing Genomic DNA was extracted from all samples at Vanderbilt University Medical Center using the Qiagen DNeasy PowerSoil HTP Kit (96-well plates) following the manufacturer's protocol, except the optional 4 °C incubations were skipped. Stool samples were thawed on ice and added directly to the kit plate. Nose, ear, and skin swabs were vortexed in tubes with 800 µL Qiagen PowerBead solution for 5 min; this PowerBead solution was then added to the kit plate. An extraction negative, which did not contain any template but was otherwise processed the same as the rest of the samples, was included on each extraction plate. To mechanically lyse the cells, plates were shaken at 20 Hz in a TissueLyser II system (Qiagen) for 20 minutes. Steps 16-33 of the kit manufacturer's protocol were performed on a QIAcube HT (Qiagen). One-step PCR targeting the V4 region of the 16S rRNA gene was performed using 515F/806R primers (56). MyTaq HS Mix (Bioline) was used to create amplicons, with the following cycling conditions: 95 °C for 2 min; 30 cycles of 95 °C for 20 sec, 50 °C for 15 sec, 72 °C for 5 min; 72 °C for 10 min; 4 °C inde nitely.
Positive PCR results were con rmed by the presence of a 400 bp band in 1% agarose gel electrophoresis; all negative controls were veri ed at this step to not have a visible band. The PCR products were cleaned and normalized using the SequalPrep Normalization Kit (Invitrogen). Samples and complementary controls (extraction negative, PCR negative, and ZymoBIOMICS Microbial Community Standard) were pooled and then cleaned using 1X AMPure XP beads. Sequencing was done on an Illumina MiSeq platform with 2 × 250 bp reads at the Vanderbilt Technologies for Applied Genomics (VANTAGE) core facility.
Alpha-and beta-diversity analyses were performed using the R package vegan (61). Prior to alpha-and beta-diversity analysis, counts were rare ed to the lowest library size, and richness, alpha-, and betaestimates were calculated. This process was repeated 400 times, and the results were averaged. Richness was estimated with the observed OTUs and Chao1 indices; alpha diversity was estimated with the Shannon and Simpson indices, which were converted into their corresponding Hill numbers (62).
Statistical testing between site alpha diversity was calculated using Mann-Whitney U or Kruskal-Wallace/Dunn's Post Hoc test where applicable. For beta-diversity analysis, counts were normalized to simple proportions, and pairwise Bray-Curtis dissimilarities were estimated. The PermANOVA (permutation-based analysis of variance) test as implemented in the Adonis function from the R package vegan was used to test for signi cance between overall microbial composition and groups of interest (i.e., S. hispidus compared to S. fulviventer and males compared to females) over 4000 permutations; results are indicated by "centroid" p-values. Homogeneity of variance within sample groups was tested using betadisper function; results are indicated by "dispersion" p-values.
Differential abundance of taxa in association with metadata categories was analyzed using DESeq2 (35). DESeq2 implements a model with raw absolute counts of each taxon with a negative binomial distribution and uses the estimated depth of sequencing of each sample to scale the (unknown) relative abundance that is the parameter of the negative binomial distribution. Compared with using either simple proportion-based normalization or rarefaction for controlling for differential sequencing depth, the DESeq2 approach provides improved sensitivity and speci city (63). Prior to DESeq2 analysis, we eliminated all taxa that were had an average number of < 10 reads, taxa with a minimum quantile mean fraction < 0.25, and taxa with a minimum quantile incidence fraction < 0.25; taxa with a normalized base mean (generated by DESeq2) < 10 were removed. Reported adjusted P values (q) values are the result of a Wald test with the Benjamini and Hochberg correction for multiple comparisons. To build alternative rankings of taxa in regard to their prevalence in one cotton rat species over the other, we also used stabsel and GeneSelector. The stabsel stability selection (64) approach aims to build the relative ranking of the predictor variables (taxa in this case) according to their importance for predicting the outcome. It does so by building multiple "base" models on random subsamples of the data. The elastic net model from the R package glmnet was used as the base feature selection method to be wrapped by the stability protocol. The ranking of taxa and their probability of being selected into the model were reported, as well as the probability cutoff corresponding to the per-family error rate that is controlled by this method. The GeneSelector package (36) was used as a stability feature ranking method that is based on a nonparametric univariate test. In brief, the same ranking method (package function RankingWilcoxon) was applied to multiple random subsamples of the full set of observations (400 replicates, sampling 50% of observations without replacement). RankingWilcoxon ranks features in each replicate according to the test statistic from Wilcoxon rank-sum test with regard to the outcome group (e.g. S. hispidus vs. S. fulviventer). Consensus ranking between replicates was then found with a Monte Carlo procedure (package function AggregateMC) and the features were reported in the order of that consensus. To account for different sequencing depth, the absolute abundance counts were normalized to simple proportions within each observation. For each feature, we also obtained several types of the effect size, such as common language effect size and rank biserial correlation. These 3 statistical analyses (DESeq2, stabsel, and GeneSelector) allowed for rigorous testing of each particular taxon of interest.

Metagenomic Library Preparation
A subset of fecal samples from 20 total male cotton rats (10 from each species), taken at days 34 and 111 within the rst cohort of cotton rats, underwent whole-metagenomic shotgun sequencing. From the same stool samples, genomic DNA was extracted using the Qiagen DNeasy PowerSoil Kit (Cat No./ID: 12888-100) by following the manufacturer's protocol (skipping the optional 4ºC incubations). In addition, a negative sample (which did not contain any template but was otherwise processed the same as the rest of the samples) and a positive control (ZymoBIOMICS Microbial Community Standard) were processed in parallel with samples and sequenced. Samples were normalized to 75 ng/ µL in 1X TE prior to library construction. Metagenomic libraries were prepared using the NEBNext® Ultra™ II FS DNA Library Prep Kit for Illumina® following the manufacturer's protocol for inputs ≤ 100 ng. Samples were fragmented at Whole Metagenomic Shotgun Sequence Analysis FastQC (65) followed by MultiQC (66) were used to examine data quality. Trimmomatic (67) was used to remove adaptors and trim low quality reads using the parameters: TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:75. Microbial communities were then pro led using MetaPhlAn2 (68). Differentially abundant bacteria were calculated using MetaPhlAn2's hclust2.py function by hierarchical clustering (based on Bray-Curtis dissimilarity) of the top 25 most abundant species according to the 90th percentile of the abundance in each clade as well as DESeq2. Functional, metabolic pro les were analyzed using Humann2, which aligns reads from UniRef (69) cluster abundances to the ChocoPhlAn (37) database.
This generates three outputs: UniRef IDs for gene families in reads per million, MetaCyc pathway coverage, and MetaCyc pathway abundances in copies per million (Supplemental File 8). To identify differential pathways between sample groups, associations between cotton rat species were identi ed by the Humann2.associate script and statistical testing using the Kruskal-Wallis H-test. Data presented (generated by Humann2.barplot script) is from pathway abundances (normalized as relative abundance) within each sample with unmapped/unintegrated pathways removed and was found statistically signi cant (p < 0.05 and q < 0.05). Superclasses distribution of identi ed MetaCyc pathways was manually generated using the online MetaCyc database.

Enumeration of Lactobacillus
Two frozen stool pellets were taken from 20 male cotton rats (10 S. hispidus, 10 S. fulviventer), weighed, and diluted to 45 mg/mL in sterile 1x PBS. Samples were rocked for 20 min on ice and resuspended manually with pipette mixing. 10 − 1 -10 − 3 serial dilutions were plated on Lactobacilli MRS agar (BD 288210) and incubated at 37ºC for 48 h. Colonies on 10 − 2 were counted, and 95 colonies grown from stool from S. hispidus and S. fulviventer each were randomly picked and inoculated into 1.2 mL MRS broth (BD 288130) in a sterile 96-deep-well plate. The plate was incubated at 37ºC for 20 h with no shaking. Cultures were gently mixed by pipetting, and 20% glycerol stocks were prepared for each culture.
Colony PCR was performed on each isolate by boiling the culture at 95ºC for 10 minutes, then using 10 µL as the template with Lactobacillus species-speci c primers (70) and MyTaq HS Red (Bioline®) with the following cycling conditions: 95 °C for 2 min; 30 cycles of 95 °C for 20 sec, 50 °C for 15 sec, 72 °C for 1 min; 72 °C for 10 min; 4 °C inde nitely. PCR reactions were spun at 3900 g for 10 min to remove any bacterial debris from the boiled template and run on 1% agarose gel to verify Lactobacillus-positive colonies. A new PCR was then repeated using the universal primers Uni331F/Uni797R (71) (following cycling conditions listed above), and puri ed PCR products were sent for Sanger sequencing. Bacterial isolate identity was determined using NCBI BLAST database.
Determination of bacterial load by qPCR DNA was extracted from an equal volume of normalized homogenates of cotton rat stool (described in Methods: Enumeration of Lactobacillus) using the DNeasy PowerSoil Kit (Qiagen). qPCR reactions were prepared in duplicate using BioRad iQ Supermix with Invitrogen Sybr Green following the manufacture's protocol. Universal eubacteria 16S rRNA primers (UniF340, UniR514) (72) equal volumes of extracted DNA, and targeted standards were used to determine copy number per gram of feces. Each qPCR plate included a corresponding extraction negative and a no-template negative control. A serial dilution of standards containing known bacterial copy numbers speci c to the primer pair were used as a standard curve as previously described. PCR reactions were run with a 15 sec 95ºC melting and 1 min 54ºC annealing step for 40 cycles. Cycle threshold (CT) values were plotted against the standard curve to determine copy number, and gures and statistical testing (unpaired T test) were generated using Prism version 8.
Isolation and culture of Lactobacillus strains from cotton rats' stool Glycerol stocks of identi ed Lactobacillus species were streaked on MRS agar plates, and a single colony was grown in culture using MRS broth. To determine growth parameters, each species was incubated at 37ºC without shaking, and growth e ciency was measured by turbidity. A growth curve was also estimated using a BioTek Synergy HTX plate reader at 37ºC for 24 h; OD 600 was measured every 10 min following a brief 3 second shake to mix culture. CFU counts were also taken during the log phase by plating a 3-fold serial dilution on MRS agar plates.

Declarations ETHICS APPROVAL AND CONSENT TO PARTICIPATE
No human subjects participated in this study. All animal procedures followed NIH and USDA guidelines and were approved by the Sigmovir Biosystems, Inc. IACUC.

AVAILABILITY OF DATA AND MATERIALS
The sequencing data will be deposited at the NCBI Short Read Archive (SRA) upon manuscript acceptance. Additionally, please contact corresponding authors for data requests.  The cotton rat microbiome as examined with 16S rRNA gene sequencing. (A) Study design includes 10 male cotton rats from both S. fulviventer and S. hispidus. Feces were taken across a 111-day period, while sample swabs of nose, ear, skin were taken at day 95. (B) Site richness (Obs. OTUs) and (C) alpha diversity (Shannon index) metrics at the OTU level indicate that the S. fulviventer microbiome was signi cantly richer and more diverse across all sites. Other alpha-diversity metrics (Chao1, Simpson) are shown in Figure S1. (D) Ordination of samples (Bray-Curtis dissimilarities, OTU level) reveals distinct microbiota compositions between feces and external sites regardless of species. Fecal microbiomes are also signi cantly different between species while external site microbiomes were not (SF=S. fulviventer, SH=S. hispidus). Statistical testing was performed using PERMANOVA between species. Each site was plotted separately in Figure S2. ns=P>0.05, *=P≤0.05, **=P≤0.01, ***=P≤0.001, ****=P≤0.0001.

Figure 2
Differential abundance of gut taxa using 16S rRNA gene sequencing. (A) Differential abundance of gut microbial taxa between S. hispidus vs. S. fulviventer that displayed signi cant differences (p < 0.05, q < 0.05, l2fc > ±0.65) between host species. The Log2FoldChange is plotted along the x-axis, with genera ranked highest in S. hispidus (black, +l2fc) to highest in S. fulviventer (grey,-l2fc) on the y-axis. Error bars represent the log2 fold change standard error; relative abundances from either S. hispidus or S. fulviventer are denoted next to the corresponding bar. (B) Probability of a gut bacterial genus being selected into a stability selection model distinguishing cotton rat species. The probability of being selected into the model is plotted along the x-axis, with top 20 ranked genera along the y-axis. (C)