Skip to main content

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

Abstract

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.

Background

The commensal microbiome can dramatically influence 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]. These include factors such as geographic regions and associated cultures and diet in humans [8], and vendor and housing facility in animal model organisms [9, 10]. In addition, recent studies have emphasized that there is a significant role of host genetics on co-evolution of the host and its associated microbiome [11,12,13]. For example, the murine genetic background is a stronger determinant of microbiome composition and structure than environmental stimuli [14]. Similarly, genetic polymorphisms, heritability, and overall host genetics in humans can also shape how commensal bacteria evolve alongside the host [15,16,17]. The microbiome has also been instrumental in predicting and protecting against severe viral disease outcomes [18, 19]. However, this burgeoning field 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) [20], influenza A virus (IAV) [21, 22], parainfluenza virus [23, 24], measles [25], human metapneumovirus [26], enterovirus [27], and human rhinovirus (HRV) [28] due to comparable human disease outcomes [29]. Cotton rats have also provided a useful model for nasal colonization studies (especially with Staphylococcus aureus) due to their human-like nasal histology [30]. 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) [31,32,33,34,35,36].

While mice have been used extensively to study viral immune responses, several factors render the mouse impractical for understanding viral pathology and kinetics, such as those relating to RSV: low replication [not translatable to humans, e.g. RSV replication is 100-fold higher in cotton rats, similar to humans [37]], resistance to upper respiratory infection [unlike humans and cotton rats, RSV does not infect the mouse nasal cavity [38, 39]], divergent lung cell infection [RSV infects ciliated bronchial epithelial cells and alveolar cells in humans and cotton rats, while only infecting pneumocytes in mice [40, 41]], and histological outcomes inconsistent with those similarly seen in the upper and lower airway of both humans and cotton rats [42]. Studies in cotton rats have also accurately predicted efficacy of several RSV therapeutics and vaccines currently used in high-risk human populations [43,44,45,46]. In light of all these factors, the cotton rat provides a superior model for studying viral-bacterial interactions than mice.

Furthermore, understanding the cotton rat microbiome is instrumental for understanding microbial interactions with viral infections, as many associational effects of both nasal and gut microbiome composition and modulation of viral outcomes have been well described in mouse models [47,48,49,50]. The microbiome of humans [9], mice [10, 51], rats [52], and other animals that are publicly available and used for answering questions relating to host microbiome and disease outcomes have been comprehensively studied and characterized. However, there has not yet been a comprehensive analysis of the microbiome in healthy cotton rats species commonly used in research (i.e., Sigmodon hispidus and S. fulviventer), making studies of viral-microbiota interactions in this animal model challenging. To date, only one study has examined the nasal microflora of healthy S. hispidus but was limited by the sample number and lack of longitudinal timepoints [53].

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, S 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 community 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 influence of host genetics on the microbiome, but also describes an appropriate animal model for studying microbiome influences 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 amplified then sequenced on the Illumina MiSeq platform with 2 × 250 base pair reads, generating an average of 35,194 reads per sample. Microbiome data was processed by following the mothur SOP, and operational taxonomic units (OTUs) were clustered at 97% identity. For diversity testing, we implemented a sample read cutoff of > 10,000/sample, which utilized 92.9% of samples for analysis with lowest library size of 11,820 reads. Remaining tests to assess differences in abundance of specific taxa (i.e., DESeq2, stabsel, GeneSelector, and LEfSE) used samples passing a per sample read cutoff of > 1000/sample, which utilized 96.4% of samples, with lowest library size of 3678 reads. We then examined the association of alpha and beta diversity (using vegan R package) and abundance of taxa (using DESeq2, stabsel, GeneSelector R packages, and LEfSE galaxy portal [54,55,56]) at each site. To compare community characteristics of cotton rats with other species, we compared beta diversity (Bray-Curtis dissimilarity) between cotton rat fecal samples from Day 0 and both humans from multiple countries [9] and mouse [51] 16S rRNA sequencing data and found that each species had a distinct community composition (Figure S1).

Fig. 1
figure1

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 significantly richer and more diverse across all sites. Other alpha-diversity metrics (Chao1, Simpson) are shown in Figure S2. d Ordination of samples (Bray-Curtis dissimilarities, OTU level) reveals distinct microbiota compositions between feces and external sites regardless of species. Color and shape of each point indicate the species and body site sampled. e-f Zoomed-in beta-diversity relationships between e fecal microbiomes and f external microbiomes. The color of each point is matched to each individual cotton rat. Statistical testing was performed using PERMANOVA between species. Each site was plotted separately in Figure S3. ns = P > 0.05, * = P ≤ 0.05, ** = P ≤ 0.01, *** = P ≤ 0.001, **** = P ≤ 0.0001

Differences in the microbiome community structure and composition between cotton rat species

Phyla within S. fulviventer and S. hispidus external sites (aggregation of ear, nose, and skin samples) were similarly dominated by Proteobacteria, Actinobacteria, and Firmicutes, with only Tenericutes being significantly more abundant in S. fulviventer (DESeq2 testing; log2 fold change = 1.04, q = 3.79E05). 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, q = 1.69E-06) and higher Firmicutes abundance in S. hispidus compared to S. fulviventer (36.2% vs. 31.8% respectively, q = 1.18E-03) (Table 1). The distribution and differential abundance of top 20 genera in each cotton rat species is shown in Figure S2.

Table 1 Percentage abundance of various phyla in both S.hispidus/S.fulviventer with standard deviations

We found that the cotton rat gut microbiome was stable over time in both cotton rat species. Richness and alpha diversity did not significantly change over time, no taxa were significantly differentially abundant when experimental day was set as the outcome variable, and beta diversity testing revealed no significant 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 (ear, nose, skin, feces) 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 S3; 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, beta dispersion p = 0.00025, Figures S4, S5D). However, comparison of beta diversity metrics of individual external sites from S. fulviventer and S. hispidus did not show significant differences (Fig. 1d; Figure S5A-C).

Fecal community

Analysis using the DESeq2 [55] package identified several bacterial genera that were differentially abundant in the two cotton rat species. These differences were most apparent in the gut (Fig. 2a). S. hispidus had a higher abundance of 18 unique genera in the gut (q < 0.05), including Lactobacillus (log2 fold change = 3.12, q = 3.27E-13), Helicobacter (log2 fold change = 2.45, q = 2.33E-33), Anaerostipes (log2 fold change = 2.35, q = 0.029), and Bifidobacterium (log2 fold change = 1.99, q = 1.03E-06). Escherichia/Shigella was more abundant in the S. hispidus gut (log2 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 (log2 fold change = − 7.19, q = 7.77E-12), Elusimicrobium (log2 fold change = − 6.22, q = 8.04E-18), and Hespellia (log2 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 confirmed using the GeneSelector [56] R package (Figure S6, Supplemental File 2). To ensure that no observed differentially abundant taxa were false-positive observations due to low abundances, we utilized the conservative LEfSe test for differential taxa [54], which reported 38 genera and confirmed 35 of 36 DESeq2-calculated differentially abundant genera (except Clostridia_unclassified; Supplemental File 3). A stability selection model showed Lactobacillus as one of the top genera (as well as those unclassified 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 S2), it was not significantly differentially abundant between the two species at any other sites except feces (Fig. 2c).

Fig. 2
figure2

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 significant 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. Unclassified phylogenetic levels indicate the lowest possible classification for that specific OTU. 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 Bacterial load depicted as copy number/uL of extracted DNA from normalized cotton rat stool. Data were generated by qPCR; and statistics were performed using unpaired T test. d Amount of aerobic colony forming units (CFU) per gram of feces on both Lactobacillus-selective (MRS) and general growth (TSA) media. S. hispidus displayed a higher amount of aerobic growth using both methods. Significance was calculated by the student’s t-test. e Percentage of CFUs with positive detection of Lactobacillus amplicons determined by PCR with primers targeting the Lactobacillus 16S rRNA region. Species-specific identity of colonies was confirmed by Sanger sequencing. “Other genera” include Bacillus, Enterococcus, and Corynebacterium. f Relative abundance of Lactobacillus between cotton rat body sites. Significance was calculated by the student’s t-test. In the gut, Lactobacillus was one of the higher-abundant taxa with significantly differential abundance between cotton rat species. Figures c-f were generated in Prism 8

External sites (nose, ear, and skin)

The following sections detail the taxa that were found to be significantly differentially abundant at each site. The full DESeq2 and GeneSelector data at both genus and family levels are shown in Supplementary Files 1 and 2 respectively. LEfSe data at the genus level is shown in Supplementary File 3. All relative abundance values at all phylogenetic levels is shown in Supplementary File 6.

Nose community

DESeq2 testing revealed that the S. hispidus nose had a higher abundance of Enterobacteriaceae (log2 fold change = 1.13, q = 1.39E-04) and Corynebacteriaceae (log2 fold change = 0.445, q = 0.00589), while the S. fulviventer nose had a higher abundance of Leuconostocaceae (log2 fold change = − 1.56, q = 0.064) (Supplementary File 1). LEfSE also revealed that S. fulviventer nose had higher abundance of Facklamia (LDA = 3.343, p = 0.0287), Bifidobacterium (LDA = 3.259, p = 0.0343), Turcibacter (LDA = 3.498, p = 0.0161), Streptococcus (LDA = 4.023, p = 0.00389), and Actinobacillus (LDA = 3.969, p = 0.00285) (Supplementary File 3).

Ear community

By DESeq2 testing, we found that the S. hispidus ear had a higher abundance of Enterobacteriaceae (log2 fold change = 0.48, q = 0.022, confirmed by LEfSe), Corynebacteriaceae (log2 fold change = 0.571, q = 0.0470), and Pseudomonas (log2 fold change = 1.33, q = 0.070, confirmed by LEfSe). DESeq2 testing showed that the S. fulviventer ear had a higher abundance of Leptotrichiaceae (log2 fold change = 0.571, q = 0.0470), Barnesiella (log2 fold change = − 4.58, q = 0.066), and Porphyromonadaceae (log2 fold change = − 2.11, q = 0.063) (Supplementary File 1). LEfSE confirmed these 3 taxa, as well as higher abundance of Sphingobacterium (LDA = 3.299, p = 0.00145), Streptococcus (LDA = 3.907, p = 0.00426), and Actinobacillus (LDA = 3.830, p = 0.00145) in S. fulviventer (Supplementary File 3).

Skin community

Enterobacteriaceae was more abundant on the S. hispidus skin (DESeq2 log2 fold change = 0.47, q = 0.0031, confirmed by LEfSe). LEfSe also found more abundant Streptococcus (LDA = 3.956, p = 0.00769), Lactococcus (LDA = 3.494, p = 0.0456), Actinobacillus (LDA = 3.730, p = 0.0129), and Mycoplasma (LDA = 4.101, p = 0.0330) on the S. fulviventer skin. (Supplementary Files 1, 3).

Confirmation 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 significantly higher in S. hispidus than S. fulviventer (Fig. 2c). Variance of bacterial load was different between the two species: while all S. fulviventer generally had a low bacterial load, the bacterial load had a large range in S. hispidus. We also plated an aliquot of normalized, homogenized stool on Lactobacillus-specific agar (De Man, Rogosa and Sharpe agar) and observed that the number of colony-forming units (CFU) per gram in S. hispidus stool was significantly higher than in S. fulviventer stool (Fig. 2d). We found that 86% of colonies picked 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 significant trend supports the relative abundance of Lactobacillus as determined by 16S rRNA gene sequencing, where Lactobacillus was significantly more abundant in S. hispidus compared to S. fulviventer (Fig. 2f). Additionally, Corynebacterium and Bacteroides species were also identified in S. hispidus stool samples. Sanger sequencing of isolates from S. fulviventer stool identified the presence of Enterococcus gallinarum and E. casselifavus.

Differences in the microbiome community structure and composition based on sex

We conducted a secondary analysis 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 significant 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 S7A-D), but there were no significant 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 S7D). Overall, differences between host sex were most pronounced in the gut compared to external sites (Figure S7) but only in S. fulviventer. Beta-diversity measurements of each species revealed that microbial composition of the gut was significantly 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 male and female 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 4.

Fig. 3
figure3

Clustering of site- and species-specific samples (Bray Curtis, OTU level) revealed host sex-dependent communities at most sites. Statistical testing was performed using PERMANOVA for both the geometric mean (or centroid) of the cluster and the dispersion (or variance). a, b Gut communities showed significant differences between both S. fulviventer (SF) and S. hispidus (SH) males and females. c-h External sites (ear, skin, and nose) showed sex-based community trends based on both sample mean distances and dispersion. Longitudinal samples from the same cotton rat are represented by matching point colors

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 first experimental group were processed for shotgun metagenomic sequencing, which generated 1.11 × 109 reads (5.57 × 107 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 classification 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 5). One sample from each S. fulviventer and S. hispidus were removed from differential expression analysis due to incongruent fitting of hierarchical clustering. Lactobacillus reuteri, L. gasseri, and the novel L. sp. ASF360 predominated the gut of S. hipsidus (Fig. 4a), and many other Lactobacillus species were significantly more abundant in S. hispidus samples compared to S. fulviventer (Figure S8). Total Lactobacillus within the S. fulviventer gut was significantly less abundant but included species unique to S. fulviventer, including L. murinus, L. BHWM-4, and L. animalis. Akkermansia muciniphilia was significantly 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 5). Proportional counts and raw counts can be found in Supplementary Files 6 and 7 respectively.

Fig. 4
figure4

Differential abundance of cotton rat gut taxa and corresponding pathways using whole-genome metagenomic sequencing. a Hierarchical clustering (Bray-Curtis) of the top 25 most differentially abundant bacterial species between S. fulviventer vs. S. hispidus. Lactobacillus reuteri and L. gasseri were drastically more abundant in S. hispidus stool, while Akkermansia muciniphila was more abundant in S. fulviventer stool. b Distribution of MetaCyc metabolomic pathways predicted from bacterial sequences. All 418 unique pathways found were represented in S. hispidus, and 334 were shared between S. hispidus and S. fulviventer. Most pathways were classified within the Biosynthesis, Superpathway, and Degradation/Utilization/Assimilation pathway superclasses. c Distribution of pathway ontology for both S. hispidus and S. fulviventer. Of all the identified pathways, the largest group consisted of Cofactor, Prosthetic Group, Electron Carrier, and Vitamin Biosynthesis, which are often components of host biology sourced solely by commensal bacteria

Differential functional potential between cotton rat species microbiome

To understand the biological implications of these differences, HUMAnN2 [57] was used to map any functional differences [MetaCyc pathway database [58]] defined by identified gene families and bacterial profiles. We identified 418 pathways (with nearly all associated with bacteria present in the sample) in the two cotton rat species based on the MetaCyc database, with all 418 pathways represented in S. hispidus but only 334 pathways represented in S. fulviventer (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 specifically, 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, p < 0.05), and most of these involved biosynthesis (Supplemental File 8).

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 identifications 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 S9). These included L-isoleucine 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 9.

Fig. 5
figure5

Several pathways that were more active in S. hispidus than S. fulviventer were greatly contributed to by Lactobacillus species. a Catalysis of proline biosynthesis by bacterial enzymes (PWY-4981). b Catalysis of the conversion of D-galactose to D-glucopyranose 6-phosphate, the more metabolically versatile carbohydrate that can feed directly into glycolysis, by the enzymes of the Leloir pathway (PWY66–422). c De novo biosynthesis of purines (PWY-6123). d De novo synthesis of ADP for the direct feeding of ATP generation, a pathway that can only accept ribonucleoside diphosphates instead of the mono- or triphosphate forms (PWY-7220). e Anaerobic breakdown of glucose to energy (PWY-5100). Table 1 Mean abundance (%) and standard deviation of the most abundant phyla identified with 16S rRNA gene sequencing in both S.fulviventer/S.hispidus (listed respectively)

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-specific 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 significantly 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 [59]. In the wild, S. hispidus and S. fulviventer are sympatric species, with S. fulviventer being the more dominant animal [60]. Separate inbreeding of the two species has made them a useful small animal model in laboratory research [29, 61]. 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. The cotton rat is a useful model for respiratory viral infection and therapeutics. In light of recent literature suggesting that the microbiome may play a key role in respiratory viral disease exacerbation or remediation and vaccine response [62,63,64,65,66], our findings exemplify the usefulness of the cotton rat model for understanding viral pathogenesis and treatments 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 [14, 67,68,69], 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 significantly higher amount of probiotic gut bacteria genera (Lactobacillus, Bifidobacterium) that have been associated with protection against the severe outcomes from RSV, IAV, and HRV [50, 70, 71]. 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 [72]. With the presence/absence of these bacteria in this animal model, along with our recently published annotated transcriptome [73], 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. In addition, our metagenomic sequencing analysis revealed species level differences as well as the metabolic potential of the microbiome. Differentially abundant taxa were directly related to differentially abundant metabolic pathways between cotton rat species. Pathways such as biosynthesis and degradation pathways (cell structures, electron carriers, vitamins, fatty acids, lipids, amino acids, etc.) could be greatly implicated in mucosal reinforcement that has been previously described in microbiome literature [74]. This microbiome-metabolic characterization may provide an important resource in understanding particular mechanisms by which the microbiome protects against certain disease states.

Our study has significant strengths compared to the lone study published so far [53], 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 first cohort of adult cotton rats. 3) We only conducted shotgun metagenomic sequencing on the gut microbiome from male cotton rats from our first cohort. While there is no cotton rat genome in the public databases, with the de novo assembly of the cotton rat transcriptome [73], further studies integrating the microbiome data and gene expression patterns may uncover more relevant information in regard 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 specifically designed for human microbiome analyses and can result in misclassification 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 platform for future hypothesis testing experiments in understanding the role of microbiome in viral pathogenesis, especially for RSV and Influenza virus. Additionally, this study adds to the small but expanding literature in understanding the role of host genetics on microbiome composition and structure.

Methods

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 first 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 fighting, 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 Ultrafine 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 min, and residual liquid was absorbed from each ear with Beaver Visitec Ultracell PVA Eye Spears pack of 5 (intended for fluid 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 mins; 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 min. 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 [75]. 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 s, 50 °C for 15 s, 72 °C for 5 min; 72 °C for 10 min; 4 °C indefinitely. Positive PCR results were confirmed by the presence of a 400 bp band in 1% agarose gel electrophoresis; all negative controls were verified 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 2x250bp reads at the Vanderbilt Technologies for Applied Genomics (VANTAGE) core facility.

16S rRNA gene data processing and statistical analysis

After sequencing, reads were processed using the mothur SOP (https://mothur.org/wiki/miseq_sop/) [76]. Operational taxonomic units (OTUs) were clustered at 97% sequence identity. Non-bacterial sequences, low-quality sequences (1.5% of total reads), and chimeras as identified with UCHIME [77] were removed during data processing. Sequences were taxonomically assigned by the Ribosomal Database Project (RDP) database 14 [78] using the SILVA database release 128 [79]. Samples with < 10,000 final reads (n = 10) were removed prior to alpha and beta diversity analysis, and samples with < 1000 final reads (n = 5) were removed prior to the remaining analyses to examine specific differentially abundant taxa (i.e., DESeq2, GeneSelector, stability selection, and LEfSE). Statistical analyses were performed using MGSAT [https://bitbucket.org/andreyto/mgsat] [18, 71], which facilitates data analysis by wrapping the R packages as described below.

Alpha- and beta-diversity analyses were performed using the R package vegan [80]. Prior to alpha- and beta-diversity analysis, counts were rarefied to the lowest library size, and richness, alpha-, and beta-estimates 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 [81]. 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 significance 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. Comparisons between Sigmodon cotton rats, human, and mouse fecal microbiome communities were performed using the same methods, and data was downloaded from NCBI Short Read Archive database (BioProject PRJNA368790, PRJEB27068, and PRJEB27068). All downloaded data was sampled from a single time point and does not represent longitudinal sampling.

Differential abundance of taxa in association with metadata categories was analyzed using DESeq2 [55] . 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 to adjust 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 [82] 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 [56] 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. LEfSe (Linear discriminant analysis [LDA] Effect Size) was executed using the online Galaxy module [54] to determine taxa most likely to explain differences between classes (species, sex, etc) using feature ranking followed by Kruskall-Wallis and pairwise Wilcoxon tests. These 4 statistical analyses (DESeq2, stabsel, GeneSelector, and LEfSe) 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 first 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 37 °C for 12 min to yield a fragment size of 200–450 bp. NEBNext Multiplex Adaptors were diluted 10-fold. NEBNext Multiplex Oligos for Illumina (Set 1, NEB #E7335) were used for PCR enrichment of adaptor-ligated DNA, and 5 cycles of PCR were run. Library quality was assessed on an Agilent 2100 Bioanalyzer System using the Agilent High Sensitivity DNA Kit (5067–4626). Samples were sequenced via the NovaSeq 6000 2 × 150 platform for Illumina at the Vanderbilt Technologies for Advanced Genomics (VANTAGE) core, aiming for 40 million reads per sample.

Whole metagenomic shotgun sequence analysis

FastQC [83] followed by MultiQC [84] were used to examine data quality. Trimmomatic [85] was used to remove adaptors and trim low quality reads using the parameters: TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:75. An average of 85% of reads mapped to various host DNA databases, but reads were not filtered before functional classification. Microbial communities were then profiled using MetaPhlAn2 [86] . 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 profiles were analyzed using HUMANn2, which aligns reads from UniRef [87] and clusters abundances to the ChocoPhlAn [57] 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 9). To identify differential pathways between sample groups, associations between cotton rat species were identified 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 significant (p < 0.05 and q < 0.05). Superclasses distribution of identified 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 were randomly picked from each species 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 min, then using 10 μL as the template with Lactobacillus species-specific primers [88] and MyTaq HS Red (Bioline®) with the following cycling conditions: 95 °C for 2 min; 30 cycles of 95 °C for 20 s, 50 °C for 15 s, 72 °C for 1 min; 72 °C for 10 min; 4 °C indefinitely. 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 [89] (following cycling conditions listed above), and purified 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) [90] 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 specific to the primer pair were used as a standard curve as previously described. PCR reactions were run with a 15 s 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 figures 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 identified 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 efficiency was measured by turbidity. A growth curve was also estimated using a BioTek Synergy HTX plate reader at 37 °C for 24 h; OD600 was measured every 10 min following a brief 3 s shake to mix culture. CFU counts were also taken during the log phase by plating a 3-fold serial dilution on MRS agar plates.

Availability of data and materials

.Sequencing data is availble through the NCBI Short Read Archive (SRA) under BioProject PRJNA721429. Additionally, please contact corresponding authors for data requests.

Abbreviations

RSV:

Respiratory Syncytial Virus

Influenza A Virus:

Influenza A Virus

HRV:

Human Rhinovirus

rRNA:

Ribosomal Ribonucleic Acid

WMS:

Whole Metagenomic Sequencing

OTU:

Operational Taxonomic Unit

LDA:

Linear Discriminant Analysis

LEFSE:

Linear Discriminant Analysis [LDA] Effect Size

CFU:

Colony-Forming Units

RT-QPCR:

Reverse Transcription Quantitative Polymerase Chain Reaction

References

  1. 1.

    Unger SA, Bogaert D. The respiratory microbiome and respiratory infections. J Inf Secur. 2017;74(Suppl 1):S84–8.

    Google Scholar 

  2. 2.

    Atherton JC, Blaser MJ. Coadaptation of helicobacter pylori and humans: ancient history, modern implications. J Clin Invest. 2009;119(9):2475–87. https://doi.org/10.1172/JCI38605.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  3. 3.

    Hooper LV, Falk PG, Gordon JI. Analyzing the molecular foundations of commensalism in the mouse intestine. Curr Opin Microbiol. 2000;3(1):79–85. https://doi.org/10.1016/S1369-5274(99)00055-7.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    David LA, Maurice CF, Carmody RN, Gootenberg DB, Button JE, Wolfe BE, et al. Diet rapidly and reproducibly alters the human gut microbiome. Nature. 2014;505(7484):559–63. https://doi.org/10.1038/nature12820.

  5. 5.

    Ley RE, Hamady M, Lozupone C, Turnbaugh PJ, Ramey RR, Bircher JS, et al. Evolution of mammals and their gut microbes. Science. 2008;320(5883):1647–51. https://doi.org/10.1126/science.1155725.

  6. 6.

    Fierer N, Hamady M, Lauber CL, Knight R. The influence of sex, handedness, and washing on the diversity of hand surface bacteria. Proc Natl Acad Sci U S A. 2008;105(46):17994–9. https://doi.org/10.1073/pnas.0807920105.

    Article  PubMed  PubMed Central  Google Scholar 

  7. 7.

    Jernberg C, Lofmark S, Edlund C, Jansson JK. Long-term ecological impacts of antibiotic administration on the human intestinal microbiota. ISME J. 2007;1(1):56–66. https://doi.org/10.1038/ismej.2007.3.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Rothschild D, Weissbrod O, Barkan E, Kurilshikov A, Korem T, Zeevi D, et al. Environment dominates over host genetics in shaping human gut microbiota. Nature. 2018;555(7695):210–5. https://doi.org/10.1038/nature25973.

  9. 9.

    Raman AS, Gehrig JL, Venkatesh S, Chang HW, Hibberd MC, Subramanian S, et al. A sparse covarying unit that describes healthy and impaired human gut microbiota development. Science. 2019;365:1–11.

  10. 10.

    Xiao L, Feng Q, Liang S, Sonne SB, Xia Z, Qiu X, et al. A catalog of the mouse gut metagenome. Nat Biotechnol. 2015;33(10):1103–8. https://doi.org/10.1038/nbt.3353.

  11. 11.

    Limborg MT, Heeb P. Special issue: coevolution of hosts and their microbiome. Genes (Basel). 2018;9(11):549. https://doi.org/10.3390/genes9110549.

  12. 12.

    Foster KR, Schluter J, Coyte KZ, Rakoff-Nahoum S. The evolution of the host microbiome as an ecosystem on a leash. Nature. 2017;548(7665):43–51. https://doi.org/10.1038/nature23292.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  13. 13.

    O'Brien PA, Webster NS, Miller DJ, Bourne DG. Host-microbe coevolution: applying evidence from model systems to complex marine invertebrate holobionts. mBio. 2019;10:1–14.

  14. 14.

    Korach-Rechtman H, Freilich S, Gerassy-Vainberg S, Buhnik-Rosenblau K, Danin-Poleg Y, Bar H, et al. Murine genetic background has a stronger impact on the composition of the gut microbiota than maternal inoculation or exposure to unlike exogenous microbiota. Appl Environ Microbiol. 2019;85(18):1–12. https://doi.org/10.1128/AEM.00826-19.

  15. 15.

    Turpin W, Espin-Garcia O, Xu W, Silverberg MS, Kevans D, Smith MI, et al. Association of host genome with intestinal microbial composition in a large healthy cohort. Nat Genet. 2016;48:1413–7.

  16. 16.

    Goodrich JK, Davenport ER, Beaumont M, Jackson MA, Knight R, Ober C, et al. Genetic determinants of the gut microbiome in UK twins. Cell Host Microbe. 2016;19(5):731–43. https://doi.org/10.1016/j.chom.2016.04.017.

  17. 17.

    Moeller AH, Caro-Quintero A, Mjungu D, Georgiev AV, Lonsdorf EV, Muller MN, et al. Cospeciation of gut microbiota with hominids. Science. 2016;353(6297):380–2. https://doi.org/10.1126/science.aaf3951.

  18. 18.

    Rosas-Salazar C, Shilts MH, Tovchigrechko A, Chappell JD, Larkin EK, Nelson KE, et al. Nasopharyngeal microbiome in respiratory syncytial virus resembles profile associated with increased childhood asthma risk. Am J Respir Crit Care Med. 2016;193(10):1180–3. https://doi.org/10.1164/rccm.201512-2350LE.

  19. 19.

    Ichinohe T, Pang IK, Kumamoto Y, Peaper DR, Ho JH, Murray TS, et al. Microbiota regulates immune defense against respiratory tract influenza a virus infection. Proc Natl Acad Sci U S A. 2011;108(13):5354–9. https://doi.org/10.1073/pnas.1019378108.

  20. 20.

    Prince GA, Hemming VG, Horswood RL, Baron PA, Chanock RM. Effectiveness of topically administered neutralizing antibodies in experimental immunotherapy of respiratory syncytial virus infection in cotton rats. J Virol. 1987;61(6):1851–4. https://doi.org/10.1128/JVI.61.6.1851-1854.1987.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Blanco JC, Pletneva LM, Wan H, Araya Y, Angel M, Oue RO, et al. Receptor characterization and susceptibility of cotton rats to avian and 2009 pandemic influenza virus strains. J Virol. 2013;87(4):2036–45. https://doi.org/10.1128/JVI.00638-12.

  22. 22.

    Ottolini MG, Blanco JCG, Eichelberger MC, Porter DD, Pletneva L, Richardson JY, et al. The cotton rat provides a useful small-animal model for the study of influenza virus pathogenesis. J Gen Virol. 2005;86(10):2823–30. https://doi.org/10.1099/vir.0.81145-0.

  23. 23.

    Ottolini MG, Porter DD, Hemming VG, Prince GA. Enhanced pulmonary pathology in cotton rats upon challenge after immunization with inactivated parainfluenza virus 3 vaccines. Viral Immunol. 2000;13(2):231–6. https://doi.org/10.1089/vim.2000.13.231.

    CAS  Article  PubMed  Google Scholar 

  24. 24.

    Ottolini MG, Porter DD, Blanco JC, Prince GA. A cotton rat model of human parainfluenza 3 laryngotracheitis: virus growth, pathology, and therapy. J Infect Dis. 2002;186(12):1713–7. https://doi.org/10.1086/345834.

    Article  PubMed  Google Scholar 

  25. 25.

    Pfeuffer J, Puschel K, Meulen V, Schneider-Schaulies J, Niewiesk S. Extent of measles virus spread and immune suppression differentiates between wild-type and vaccine strains in the cotton rat model (Sigmodon hispidus). J Virol. 2003;77(1):150–8. https://doi.org/10.1128/JVI.77.1.150-158.2003.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Hamelin ME, Yim K, Kuhn KH, Cragin RP, Boukhvalova M, Blanco JC, et al. Pathogenesis of human metapneumovirus lung infection in BALB/c mice and cotton rats. J Virol. 2005;79(14):8894–903. https://doi.org/10.1128/JVI.79.14.8894-8903.2005.

  27. 27.

    Patel MC, Wang W, Pletneva LM, Rajagopala SV, Tan Y, Hartert TV, et al. Enterovirus D-68 infection, prophylaxis, and vaccination in a novel permissive animal model, the cotton rat (Sigmodon hispidus). PLoS One. 2016;11(11):e0166336. https://doi.org/10.1371/journal.pone.0166336.

  28. 28.

    Blanco JC, Core S, Pletneva LM, March TH, Boukhvalova MS, Kajon AE. Prophylactic antibody treatment and intramuscular immunization reduce infectious human rhinovirus 16 load in the lower respiratory tract of challenged cotton rats. Trials Vaccinol. 2014;3:52–60. https://doi.org/10.1016/j.trivac.2014.02.003.

    Article  PubMed  PubMed Central  Google Scholar 

  29. 29.

    Boukhvalova MS, Prince GA, Blanco JC. The cotton rat model of respiratory viral infections. Biologicals. 2009;37(3):152–9. https://doi.org/10.1016/j.biologicals.2009.02.017.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  30. 30.

    Burian M, Rautenberg M, Kohler T, Fritz M, Krismer B, Unger C, et al. Temporal expression of adhesion factors and activity of global regulators during establishment of Staphylococcus aureus nasal colonization. J Infect Dis. 2010;201(9):1414–21. https://doi.org/10.1086/651619.

  31. 31.

    Carrara AS, Coffey LL, Aguilar PV, Moncayo AC, Da Rosa AP, Nunes MR, et al. Venezuelan equine encephalitis virus infection of cotton rats. Emerg Infect Dis. 2007;13(8):1158–65. https://doi.org/10.3201/eid1308.061157.

  32. 32.

    Rollin PE, Ksiazek TG, Elliott LH, Ravkov EV, Martin ML, Morzunov S, et al. Isolation of black creek canal virus, a new hantavirus from Sigmodon hispidus in Florida. J Med Virol. 1995;46(1):35–9. https://doi.org/10.1002/jmv.1890460108.

  33. 33.

    Holsomback TS, McIntyre NE, Nisbett RA, Strauss RE, Chu YK, Abuzeineh AA, et al. Bayou virus detected in non-oryzomyine rodent hosts: an assessment of habitat composition, reservoir community structure, and marsh rice rat social dynamics. J Vector Ecol. 2009;34(1):9–21. https://doi.org/10.1111/j.1948-7134.2009.00003.x.

  34. 34.

    Winn WC Jr, Murphy FA. Tamiami virus infection in mice and cotton rats. Bull World Health Organ. 1975;52(4-6):501–6.

    PubMed  PubMed Central  Google Scholar 

  35. 35.

    Dietrich G, Montenieri JA, Panella NA, Langevin S, Lasater SE, Klenk K, et al. Serologic evidence of west nile virus infection in free-ranging mammals, Slidell, Louisiana, 2002. Vector Borne Zoonotic Dis. 2005;5(3):288–92. https://doi.org/10.1089/vbz.2005.5.288.

  36. 36.

    Ehlen L, Todtmann J, Specht S, Kallies R, Papies J, Muller MA, et al. Epithelial cell lines of the cotton rat (Sigmodon hispidus) are highly susceptible in vitro models to zoonotic Bunya-, Rhabdo-, and Flaviviruses. Virol J. 2016;13(1):74. https://doi.org/10.1186/s12985-016-0531-5.

  37. 37.

    Prince GA, Jenson AB, Horswood RL, Camargo E, Chanock RM. The pathogenesis of respiratory syncytial virus infection in cotton rats. Am J Pathol. 1978;93(3):771–91.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. 38.

    Gitiban N, Jurcisek JA, Harris RH, Mertz SE, Durbin RK, Bakaletz LO, et al. Chinchilla and murine models of upper respiratory tract infections with respiratory syncytial virus. J Virol. 2005;79(10):6035–42. https://doi.org/10.1128/JVI.79.10.6035-6042.2005.

  39. 39.

    Graham BS, Perkins MD, Wright PF, Karzon DT. Primary respiratory syncytial virus infection in mice. J Med Virol. 1988;26(2):153–62. https://doi.org/10.1002/jmv.1890260207.

    CAS  Article  PubMed  Google Scholar 

  40. 40.

    Johnson JE, Gonzales RA, Olson SJ, Wright PF, Graham BS. The histopathology of fatal untreated human respiratory syncytial virus infection. Mod Pathol. 2007;20(1):108–19. https://doi.org/10.1038/modpathol.3800725.

    CAS  Article  PubMed  Google Scholar 

  41. 41.

    Martinez-Sobrido L, Gitiban N, Fernandez-Sesma A, Cros J, Mertz SE, Jewell NA, et al. Protection against respiratory syncytial virus by a recombinant Newcastle disease virus vector. J Virol. 2006;80(3):1130–9. https://doi.org/10.1128/JVI.80.3.1130-1139.2006.

  42. 42.

    Grieves JL, Yin Z, Durbin RK, Durbin JE. Acute and chronic airway disease after human respiratory syncytial virus infection in cotton rats (Sigmodon hispidus). Comp Med. 2015;65(4):315–26.

    CAS  PubMed  PubMed Central  Google Scholar 

  43. 43.

    Rodriguez WJ, Gruber WC, Welliver RC, Groothuis JR, Simoes EA, Meissner HC, et al. Respiratory syncytial virus (RSV) immune globulin intravenous therapy for RSV lower respiratory tract infection in infants and young children at high risk for severe RSV infections: respiratory syncytial virus immune globulin study group. Pediatrics. 1997;99(3):454–61. https://doi.org/10.1542/peds.99.3.454.

  44. 44.

    Prince GA, Curtis SJ, Yim KC, Porter DD. Vaccine-enhanced respiratory syncytial virus disease in cotton rats following immunization with lot 100 or a newly prepared reference vaccine. J Gen Virol. 2001;82(12):2881–8. https://doi.org/10.1099/0022-1317-82-12-2881.

    CAS  Article  PubMed  Google Scholar 

  45. 45.

    Prince GA, Jenson AB, Hemming VG, Murphy BR, Walsh EE, Horswood RL, et al. Enhancement of respiratory syncytial virus pulmonary pathology in cotton rats by prior intramuscular inoculation of formalin-inactiva ted virus. J Virol. 1986;57(3):721–8. https://doi.org/10.1128/JVI.57.3.721-728.1986.

  46. 46.

    Ottolini MG, Curtis SR, Mathews A, Ottolini SR, Prince GA. Palivizumab is highly effective in suppressing respiratory syncytial virus in an immunosuppressed animal model. Bone Marrow Transplant. 2002;29(2):117–20. https://doi.org/10.1038/sj.bmt.1703326.

    CAS  Article  PubMed  Google Scholar 

  47. 47.

    Chiba E, Tomosada Y, Vizoso-Pinto MG, Salva S, Takahashi T, Tsukida K, et al. Immunobiotic Lactobacillus rhamnosus improves resistance of infant mice against respiratory syncytial virus infection. Int Immunopharmacol. 2013;17(2):373–82. https://doi.org/10.1016/j.intimp.2013.06.024.

  48. 48.

    Fonseca W, Lucey K, Jang S, Fujimura KE, Rasky A, Ting HA, et al. Lactobacillus johnsonii supplementation attenuates respiratory viral infection via metabolic reprogramming and immune cell modulation. Mucosal Immunol. 2017;10(6):1569–80. https://doi.org/10.1038/mi.2017.13.

  49. 49.

    Hyde ER, Petrosino JF, Piedra PA, Camargo CA Jr, Espinola JA, Mansbach JM. Nasopharyngeal Proteobacteria are associated with viral etiology and acute wheezing in children with severe bronchiolitis. J Allergy Clin Immunol. 2014;133(4):1220–2. https://doi.org/10.1016/j.jaci.2013.10.049.

    Article  PubMed  Google Scholar 

  50. 50.

    Tomosada Y, Chiba E, Zelaya H, Takahashi T, Tsukida K, Kitazawa H, et al. Nasally administered Lactobacillus rhamnosus strains differentially modulate respiratory antiviral immune responses and induce protection against respiratory syncytial virus infection. BMC Immunol. 2013;14(1):40. https://doi.org/10.1186/1471-2172-14-40.

  51. 51.

    Casero D, Gill K, Sridharan V, Koturbash I, Nelson G, Hauer-Jensen M, et al. Space-type radiation induces multimodal responses in the mouse gut microbiome and metabolome. Microbiome. 2017;5(1):105. https://doi.org/10.1186/s40168-017-0325-z.

  52. 52.

    Piccolo BD, Graham JL, Stanhope KL, Nookaew I, Mercer KE, Chintapalli SV, et al. Diabetes-associated alterations in the cecal microbiome and metabolome are independent of diet or environment in the UC Davis type 2 diabetes mellitus rat model. Am J Physiol Endocrinol Metab. 2018;315(5):E961–72. https://doi.org/10.1152/ajpendo.00203.2018.

  53. 53.

    Chaves-Moreno D, Plumeier I, Kahl S, Krismer B, Peschel A, Oxley AP, et al. The microbial community structure of the cotton rat nose. Environ Microbiol Rep. 2015;7(6):929–35. https://doi.org/10.1111/1758-2229.12334.

  54. 54.

    Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, et al. Metagenomic biomarker discovery and explanation. Genome Biol. 2011;12(6):R60. https://doi.org/10.1186/gb-2011-12-6-r60.

  55. 55.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. https://doi.org/10.1186/s13059-014-0550-8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  56. 56.

    Boulesteix AL, Slawski M. Stability and aggregation of ranked gene lists. Brief Bioinform. 2009;10(5):556–68. https://doi.org/10.1093/bib/bbp034.

    CAS  Article  PubMed  Google Scholar 

  57. 57.

    Franzosa EA, McIver LJ, Rahnavard G, Thompson LR, Schirmer M, Weingart G, et al. Species-level functional profiling of metagenomes and metatranscriptomes. Nat Methods. 2018;15(11):962–8. https://doi.org/10.1038/s41592-018-0176-y.

  58. 58.

    Caspi R, Billington R, Fulcher CA, Keseler IM, Kothari A, Krummenacker M, et al. The MetaCyc database of metabolic pathways and enzymes. Nucleic Acids Res. 2018;46(D1):D633–9. https://doi.org/10.1093/nar/gkx935.

  59. 59.

    Henson DD, Bradley RD. Molecular systematics of the genus Sigmodon: results from mitochondrial and nuclear gene sequences. Can J Zool. 2009;87(3):211–20. https://doi.org/10.1139/Z09-005.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  60. 60.

    Petersen MK. Interactions between the cotton rats, sigmodon fulviventer and S. hispidus. Am Midland Nat. 1973;90(2):319–33.

    Article  Google Scholar 

  61. 61.

    Sadowski W, Semkow R, Wilczynski J, Krus S, Kantoch M. [The cotton rat (Sigmodon hispidus) as an experimental model for studying viruses in human respiratory tract infections. I. Para-influenza virus type 1, 2 and 3, adenovirus type 5 and RS virus]. Med Dosw Mikrobiol. 1987;39:33–42.

    CAS  PubMed  Google Scholar 

  62. 62.

    McAleer JP, Kolls JK. Contributions of the intestinal microbiome in lung immunity. Eur J Immunol. 2018;48(1):39–49. https://doi.org/10.1002/eji.201646721.

    CAS  Article  PubMed  Google Scholar 

  63. 63.

    Ding T, Song T, Zhou B, Geber A, Ma Y, Zhang L, et al. Microbial composition of the human nasopharynx varies according to influenza virus type and vaccination status. mBio. 2019;10:1–15.

  64. 64.

    Ciabattini A, Olivieri R, Lazzeri E, Medaglini D. Role of the microbiota in the modulation of vaccine immune responses. Front Microbiol. 2019;10:1305. https://doi.org/10.3389/fmicb.2019.01305.

    Article  PubMed  PubMed Central  Google Scholar 

  65. 65.

    Hagan T, Cortese M, Rouphael N, Boudreau C, Linde C, Maddur MS, et al. Antibiotics-driven gut microbiome perturbation alters immunity to vaccines in humans. Cell. 2019;178:1313–1328.e13.

  66. 66.

    Harris VC, Armah G, Fuentes S, Korpela KE, Parashar U, Victor JC, et al. Significant correlation between the infant gut microbiome and rotavirus vaccine response in rural Ghana. J Infect Dis. 2017;215(1):34–41. https://doi.org/10.1093/infdis/jiw518.

  67. 67.

    Hufeldt MR, Nielsen DS, Vogensen FK, Midtvedt T, Hansen AK. Variation in the gut microbiota of laboratory mice is related to both genetic and environmental factors. Comp Med. 2010;60(5):336–47.

    CAS  PubMed  PubMed Central  Google Scholar 

  68. 68.

    Org E, Parks BW, Joo JW, Emert B, Schwartzman W, Kang EY, et al. Genetic and environmental control of host-gut microbiota interactions. Genome Res. 2015;25(10):1558–69. https://doi.org/10.1101/gr.194118.115.

  69. 69.

    Tabrett A, Horton MW. The influence of host genetics on the microbiome. F1000Res. 2020;9:1–9.

  70. 70.

    Kumpu M, Kekkonen RA, Korpela R, Tynkkynen S, Jarvenpaa S, Kautiainen H, et al. Effect of live and inactivated Lactobacillus rhamnosus GG on experimentally induced rhinovirus colds: randomised, double blind, placebo-controlled pilot trial. Benef Microbes. 2015;6(5):631–9. https://doi.org/10.3920/BM2014.0164.

  71. 71.

    Rosas-Salazar C, Shilts MH, Tovchigrechko A, Schobel S, Chappell JD, Larkin EK, et al. Nasopharyngeal Lactobacillus is associated with a reduced risk of childhood wheezing illnesses following acute respiratory syncytial virus infection in infancy. J Allergy Clin Immunol. 2018;142:1447–1456.e9.

  72. 72.

    Kuss SK, Best GT, Etheredge CA, Pruijssers AJ, Frierson JM, Hooper LV, et al. Intestinal microbiota promote enteric virus replication and systemic pathogenesis. Science. 2011;334(6053):249–52. https://doi.org/10.1126/science.1211057.

  73. 73.

    Rajagopala SV, Singh H, Patel MC, Wang W, Tan Y, Shilts MH, et al. Cotton rat lung transcriptome reveals host immune response to respiratory syncytial virus infection. Sci Rep. 2018;8(1):11318. https://doi.org/10.1038/s41598-018-29374-x.

  74. 74.

    Hiippala K, Jouhten H, Ronkainen A, Hartikainen A, Kainulainen V, Jalanka J, et al. The potential of gut commensals in reinforcing intestinal barrier function and alleviating inflammation. Nutrients. 2018;10(8):1–23. https://doi.org/10.3390/nu10080988.

  75. 75.

    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. https://doi.org/10.1128/AEM.01043-13.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  76. 76.

    Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, 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. https://doi.org/10.1128/AEM.01541-09.

  77. 77.

    Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. UCHIME improves sensitivity and speed of chimera detection. Bioinformatics. 2011;27(16):2194–200. https://doi.org/10.1093/bioinformatics/btr381.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  78. 78.

    Cole JR, Wang Q, Cardenas E, Fish J, Chai B, Farris RJ, et al. The ribosomal database project: improved alignments and new tools for rRNA analysis. Nucleic Acids Res. 2009;37(Database):D141–5. https://doi.org/10.1093/nar/gkn879.

  79. 79.

    Pruesse E, Quast C, Knittel K, Fuchs BM, Ludwig W, Peplies J, et al. 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. https://doi.org/10.1093/nar/gkm864.

  80. 80.

    ari Oksanen FGB, Michael Friendly, Roeland Kindt, Pierre Legendre, Dan McGlinn, Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos, M. Henry H. Stevens, Eduard Szoecs and Helene Wagner. 2019. vegan: community ecology package. R package version 2.5–6. https://CRAN.R-project.org/package=vegan. Accessed 24 Sept 2020.

  81. 81.

    Hill MO. Diversity and evenness: a unifying notation and its consequences: Ecology; Ecological Society of America; 1973.

    Google Scholar 

  82. 82.

    Nicolai Meinshausen PB. Stability selection. J R Stat Soc. 2010;72(4):417–73. https://doi.org/10.1111/j.1467-9868.2010.00740.x.

    Article  Google Scholar 

  83. 83.

    Andrews S. 2015. FastQC: a quality control tool for high throughput sequence data [online]. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/. Accessed 24 Sept 2020.

    Google Scholar 

  84. 84.

    Ewels P, Magnusson M, Lundin S, Kaller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;32(19):3047–8. https://doi.org/10.1093/bioinformatics/btw354.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  85. 85.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://doi.org/10.1093/bioinformatics/btu170.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  86. 86.

    Segata N, Waldron L, Ballarini A, Narasimhan V, Jousson O, Huttenhower C. Metagenomic microbial community profiling using unique clade-specific marker genes. Nat Methods. 2012;9(8):811–4. https://doi.org/10.1038/nmeth.2066.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  87. 87.

    Suzek BE, Wang Y, Huang H, McGarvey PB, Wu CH, UniProt C. UniRef clusters: a comprehensive and scalable alternative for improving sequence similarity searches. Bioinformatics. 2015;31(6):926–32. https://doi.org/10.1093/bioinformatics/btu739.

    CAS  Article  PubMed  Google Scholar 

  88. 88.

    Nadkarni MA, Martin FE, Jacques NA, Hunter N. Determination of bacterial load by real-time PCR using a broad-range (universal) probe and primers set. Microbiology (Reading). 2002;148(1):257–66. https://doi.org/10.1099/00221287-148-1-257.

    CAS  Article  Google Scholar 

  89. 89.

    Delroisse JM, Boulvin AL, Parmentier I, Dauphin RD, Vandenbol M, Portetelle D. Quantification of Bifidobacterium spp. and Lactobacillus spp. in rat fecal samples by real-time PCR. Microbiol Res. 2008;163(6):663–70. https://doi.org/10.1016/j.micres.2006.09.004.

    CAS  Article  PubMed  Google Scholar 

  90. 90.

    Barman M, Unold D, Shifley K, Amir E, Hung K, Bos N, et al. Enteric salmonellosis disrupts the microbial ecology of the murine gastrointestinal tract. Infect Immun. 2008;76(3):907–15. https://doi.org/10.1128/IAI.01432-07.

Download references

Acknowledgements

We thank Dr. Stokes Peebles Jr. for careful reading of the manuscript.

Funding

This work was supported by start-up funds from Vanderbilt University Medical Center awarded to SRD and funds from the National Institute of Allergy and Infectious Diseases (under award numbers 1R21AI149262 and 1R21AI154016). This work was also supported by funds from National Institutes of Allergy and Infectious Diseases to JCGB (under the award No R01AI125215). SRD is also supported by 1R21AI142321, R21AI142321-01A1S1, U19AI095227, and U19AI110819), and the National Heart, Lung and Blood Institute (1R01HL146401) the Edward P. Evans Foundation, the Vanderbilt Institute for Clinical and Translational Research (grant support from the National Center for Advancing Translational Sciences under award number UL1TR000445), Vanderbilt Technologies for Advanced Genomics Core (grant support from the National Institutes of Health under award numbers UL1RR024975, P30CA68485, P30EY08126, and G20RR030956) and Sigmovir Biosystems Inc’s corporate funds.

Author information

Affiliations

Authors

Contributions

SRD, JCGB, BAS conceived the idea. MCP, AK, WZ, DS, and MSB maintained cotton rat colonies and collected samples under the supervision of JCGB. BAS and HHB performed DNA extraction and 16S rRNA gene sequencing library preparation. BAS performed whole metagenomic sequencing library preparation, bacterial cultivation, and qPCR. BAS, MHS, SVR performed sequencing data analysis. MHS, SY, SVR, SRD were consulted during data analysis. CRS assisted in manuscript preparation and data presentation oversight. BAS wrote the manuscript with inputs from SRD, MHS and CRS. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Jorge C. G. Blanco or Suman R. Das.

Ethics 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.

Consent for publication

Not Applicable.

Competing interests

JCGB, MCP, AK, WZ, DS, and MSB are employed by Sigmovir Biosystems, Inc. Other authors have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Ordination of human, mouse, and two Sigmodon cotton rat species (Bray-Curtis dissimilarities, OTU level) reveals that the cotton rat fecal microbiome is very distinct from that of humans, and more similar to, but still distinct from, mice.

Additional file 2: Figure S2.

Top 20 most abundant bacterial genera at each body site of S. hispidus and S. fulviventer. In both species of cotton rats, external sites (skin, ear, nose) shared similar dominating genera while there were notable difference in gut taxa between S. hispidus and S. fulviventer. Not all reads were able to be classified down to the genus level; the lowest taxonomic level available is reported. The letter after the classification denotes the lowest taxonomic level able to be identified for the particular OTU (i.e., g for genus, f for family, o for order, c for class, p for phylum).

Additional file 3: Figure S3.

Chao1 and Simpson alpha diversity indices of ear, nose, skin, and feces between S. fulviventer and S. hispidus. Statistical testing between each cotton rat species was performed using a Student’s t-test. Statistical testing between body sites is not shown (no significant differences across external sites). ns = P > 0.05, * = P ≤ 0.05, ** = P ≤ 0.01, *** = P ≤ 0.001, **** = P ≤ 0.0001.

Additional file 4: Figure S4.

PCoA plots comparing Bray-Curtis dissimilarities between body sites in both (A) S. hispidus and (B) S. fulviventer.

Additional file 5: Figure S5.

Difference in body site beta diversity between S. hispidus and S. fulviventer. Clustering of samples (Bray-Curtis, OTU level) shows separation by host species (A) Ear, (B) Skin, (C) Nose, (D) Feces.

Additional file 6: Figure S6.

Statistically significant differentially abundant bacteria genera between S. hispidus and S. fulviventer determined by GeneSelector.

Additional file 7: Figure S7.

Alpha diversity metrics of ear, nose, skin, and feces between male and female S. fulviventer and S. hispidus. Richness and diversity were determined using the following methods: A) Observed OTUs, B) Chao1 index, C) Shannon Diversity Index, and D) Simpson Diversity Index. Statistical testing between gender of each cotton rat species was performed using a one-way (feces) and two-way (external sites) ANOVA, followed by a Tukey’s post-hoc test. The p-value as a result of all comparisons is shown in the top right; pairwise comparisons that were found to be significant with the Tukey post-hoc test are denoted by a bar and asterisk above the groups being compared. Statistical testing across species is not shown. Feces were plotted separately to account for the discrepancy between the Y axes. ns = P > 0.05, * = P ≤ 0.05, ** = P ≤ 0.01, *** = P ≤ 0.001, **** = P ≤ 0.0001.

Additional file 8: Figure S8.

DESeq2 results representing differential abundance of Lactobacillus species and strains between S. hispidus and S. fulviventer. Data were generated from whole metagenome sequencing.

Additional file 9: Figure S9.

Several pathways that were more active in S. fulviventer than S. hispidus were greatly contributed to by Akkermansia species. (A) A member of the superpathway of branched chain amino acid biosynthesis, that generates not only isoleucine, but also leucine and valine (ILEUSYN-PWY). (B) Degradation of starch for the generation of carbon skeletons, reductants, and ATP for anabolic bacterial fatty acid pathway initiation via pyruvate decarboxylation to acetyl CoA (PWY-1042). (C) Biosynthesis of R-4′-phosphopantothenate, the universal precursor for the synthesis of coenzyme A and acyl carrier protein. Only plants and microorganisms can synthesize pantothenate de novo; animals require a dietary supplement. Synonymous with Vitamin B5 synthesis. (PANTO-PWY). (D) A member of the superpathway of branched chain amino acid biosynthesis, that generates not only valine, but also leucine and isoleucine.

Additional file 10: Table S1.

Differential abundance analysis of taxa between individual body site across female and male S. hispidus and S. fulviventer. Positive Log2FoldChange = higher in females; negative Log2FoldChange higher in males. There were no significant taxa in S. fulviventer feces.

Additional file 11: Supplemental File 1.

DESeq2 analysis of 16S rRNA gene sequencing data between S. fulviventer and S. hispidus across all body sites at the genus and family level.

Additional file 12: Supplemental File 2.

GeneSelector analysis of 16S rRNA gene sequencing data between S. fulviventer and S. hispidus across all body sites at the genus and family level.

Additional file 13: Supplemental File 3.

LEfSe analysis of 16S rRNA gene sequencing data between S. fulviventer and S. hispidus across all body sites at the genus level.

Additional file 14: Supplemental File 4.

DESeq2 analysis of 16S rRNA gene sequencing data between male and female cotton rats (both S. fulviventer and S. hispidus) across all body sites at the genus and family level.

Additional file 15: Supplemental File 5.

DESeq2 analysis of whole metagenomic sequencing data between S. fulviventer and S. hispidus across all body sites at the species level.

Additional file 16: Supplemental File 6.

Proportional counts from 16S rRNA sequencing data at all levels for both cotton rat cohorts.

Additional file 17: Supplemental File 7.

Raw counts from 16S rRNA sequencing data at all levels for both cotton rat cohorts.

Additional file 18: Supplemental File 8.

All annotated MetaCyc pathways identified by HUMANn2 from whole metagenomic sequencing data.

Additional file 19: Supplemental File 9.

Statistical comparison of MetaCyc pathway comparison between S. fulviventer and S. hispidus.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Strickland, B.A., Patel, M.C., Shilts, M.H. et al. Microbial community structure and composition is associated with host species and sex in Sigmodon cotton rats. anim microbiome 3, 29 (2021). https://doi.org/10.1186/s42523-021-00090-8

Download citation

Keywords

  • Microbiome
  • Metagenomics
  • Cotton rat
  • Sigmodon
  • 16S rRNA gene
  • Gut
  • Skin
  • Respiratory
\