Biodiversity indices and Random Forests reveal the potential for striped skunk (Mephitis mephitis) fecal microbial communities to function as a biomarker for oral rabies vaccination

Abstract

Wildlife disease surveillance and monitoring poses unique challenges when assessing rates of population vaccination, immunity, or infection prevalence. Non-invasively detected biomarkers can help reduce risk to both animal and field personnel during wildlife disease management activities. In this study, we investigated the utility of fecal microbiome data collected from captive striped skunks (Mephitis mephitis) in predicting rabies virus vaccination and infection status. We sequenced the hypervariable region 4 (V4) of the bacterial 16S gene and estimated alpha and beta diversity across timepoints in three groups of skunks: vaccination then rabies virus infection, sham vaccination then rabies virus infection, and rabies virus infected without vaccination. Alpha diversity did not differ among treatment groups but beta diversity between treatments was statistically significant. The phyla Firmicutes and Proteobacteria were dominant among all samples. Using Random Forests, we identified operational taxonomic units (OTUs) that greatly influenced classification of fecal samples into treatment groups. Each of these OTUs was correlated with fecal volatile organic compounds detected from the samples for companion treatment groups in another study. This research is the first to highlight striped skunk microbiome biodiversity as a vaccination biomarker which pushes the frontier on alternative methods for surveillance and monitoring of vaccination and disease in wildlife populations.

Introduction

Most of our knowledge of non-human animal microbiomes comes from studies on laboratory animals or livestock, with a handful of studies on wildlife [1, 2]. However, there is now greater recognition of the role wildlife-specific microbiomes play in understanding and managing wildlife and associated pathogens [3, 4]. The microbiome is a crucial part of the mammalian immune system which can help determine infection outcomes by preventing the establishment of pathogens or facilitating immune responses [5]. Pathogens and other immune system triggers, such as vaccination, can induce concomitant responses in microbiome community structure. In some cases, immune system activation has predictable outcomes on the microbiome which has led some to suggest that alterations to taxonomic diversity can be used as a host biomarker of pathogen exposure, vaccination, and immune response [e.g., 6–9].

In wildlife disease management, one challenge is determining the extent of infection, vaccination, and immunity in natural populations. Biomarkers can inform monitoring and surveillance by providing vaccination estimates and non-invasive biomarkers could greatly ease the burden on agencies involved with animal trapping for disease surveillance and management. Biomarkers can include volatile organic compounds (VOCs) excreted by wildlife in urine, feces, and breath, and have been used to detect pathogen infections or vaccination status [e.g., 10–13]. Given that VOCs are often byproducts of microbial physiology, one can expect that changes in microbial communities may be correlated with changes in VOC profiles. Thus, if VOCs can be utilized as biomarkers of pathogen infection and vaccination, then the underlying host microbiome may also serve this purpose.

Shifts in microbial communities of wild mammals in response to environmental influences have been demonstrated for stressors such as habitat fragmentation, pollution, disease, and hibernation [3, 4, 14, 15]. While researchers have been able to detect shifts in microbiome communities, one question often remains: can certain microbial communities predict the presence or exposure to a immunological stressor, such as a pathogen or vaccine? Abdo, LeCureux [7] identified bacterial operational taxonomic units (OTUs) in lab mice that were associated with treatment groups in a study on a probiotic oral vaccine vector for human immunodeficiency virus (HIV). In two similar studies, Zhang, Wu [6] also found bacterial taxa that were associated with rabies virus infection and oral rabies vaccination in lab mice, and Shi, Zou [16] found segmented filamentous bacteria (Candidatus arthromitus) helps prevent rotavirus infection and alleviates infection-related symptoms in a murine model.

While little is known about the use of microbiomes as biomarkers, they hold great potential for utility in wildlife disease surveillance. One pertinent example of using the microbiome to detect pathogen infection was the identification of cloacal bacterial OTUs that were associated with influenza A virus infection status in wild birds [17]. Collecting fecal microbiome samples from wildlife is considered a non-invasive sampling technique that reduces the burden and risk to the animal and field personnel while allowing for collection of robust sample sizes [18]. Further, individual genetic identification can be obtained from fecal samples, and when used in conjunction with mark-recapture approaches, the population proportion of animals that are infected or vaccinated could be estimated.

Rabies lyssavirus (i.e., rabies virus) is a highly lethal zoonotic pathogen that can infect all mammal species. Rabies virus predominantly occurs in Chiropteran and Carnivoran wildlife reservoirs, with co-evolved variants that circulate within a single species with occasional spillover events to other mammals [19–21]. On the North American continent, rabies virus lineages of bat-origin and domestic dog-origin (Canis lupus familiaris) circulate in wild mesocarnivore species [21]. Within these clades there are subclades of variants that are found primarily in specific species [22, 23]. In the United States (U.S.), the United States Department of Agriculture (USDA), along with other Federal, state, and local agencies, work to mitigate the risks posed to humans, domestic animals, livestock and wildlife from rabies virus, by preventing the spread and local elimination of the raccoon (Procyon lotor) variant of rabies virus in the eastern U.S., and maintaining a canine rabies virus free status, through a program of oral vaccination and enhanced surveillance [24, 25]. In the central and western U.S., striped skunks (Mephitis mephitis) are one of the principal wildlife reservoirs for rabies virus consisting of three variants: northcentral skunk, southcentral skunk, and California skunk [23]. Striped skunks have also been recognized as the most frequent spillover recipient of raccoon variant rabies virus and are a key secondary target of control efforts in the eastern U.S. [26, 27].

The USDA currently conducts enhanced surveillance of rabies virus in mesocarnivores to track changes in epizootiology and to evaluate the effectiveness of vaccination programs for control and elimination [28]. Sampling requires trapping and handling of animals, which are potentially rabid. Field staff may also check for oral vaccine bait consumption by removing the premolar of captured raccoons; where vaccine-baits contain tetracycline which incorporates into the teeth functioning as a biomarker [29]. A non-invasive surveillance and monitoring tool could be a useful alternative to current strategies by alleviating the need to capture or handle animals to remove brain tissues or teeth respectively. Two studies have demonstrated that fecal samples may be an effective way to detect rabies virus infection and vaccination. Kimball, Volker [12] determined that certain fecal volatile organic compounds (VOCs) collected from captive striped skunks were correlated with rabies virus binding antibody response to oral rabies vaccination. Zhang, Wu [6] demonstrated that rabies virus neutralizing antibody (RVNA) response to vaccination was degraded following antibiotic treatment of lab mice and that certain bacterial taxa detected in fecal samples were positively or negatively correlated with high and low RVNA titers. In this follow up study to Kimball, Volker [12] we collected fecal samples from captive striped skunks used in oral rabies vaccine efficacy trials, to evaluate whether vaccination and/or rabies virus infection alters the fecal microbiome and thus could be used as an indicator of vaccination status.

Methods

Ethics statement

Animal use and procedures described were consistent with the guidelines of the USDA National Wildlife Research Center (NWRC) Institutional Animal Care and Use Committee (protocol QA-2258) and import and housing of skunks were permitted under Colorado Parks and Wildlife 13TR2056A1, 14TR2056, and 15TR2143.

Vaccine trials and fecal collection

The animal handling and experimental treatments are described in detail in Gilbert, Johnson [30] and Kimball, Volker [12] but we provide a brief summary of vaccination, rabies virus infection, and fecal collection. We included samples from 22 animals in this study which was a subset of the 35 in Gilbert, Johnson [30] and included all animals from Kimball, Volker [12]. The vaccinated individuals were randomly assigned to one of four groups. Three of the groups received three different doses of Ontario Rabies Vaccine, ONRAB (Artemis Technologies, Inc., Guelph, Ontario, Canada) diluted in minimal essential media (MEM) supplemented with 5% fetal bovine serum (10^10.2, 10^9.8, or 10^9.3 median tissue culture infective doses (TCID₅₀; all combined into group Rvax for this study)). Ten skunks were given the oral rabies vaccine (group Rvax), five animals were given sham vaccine orally (MEM only; Group SHvax), and seven skunks were unvaccinated and challenged with rabies virus (no vaccine; group Rinf). The vaccine or sham treatments were administered to individual skunks by direct instillation into the oral cavity under light anesthesia with isoflurane gas [30]. On day 335 post-vaccination skunks were challenged with a New York City dog variant of rabies virus (92-5A) as described in Gilbert, Johnson [30].

Control skunks (group Rinf) were inoculated with 92-5A rabies virus as described in [30]. These individuals were randomly sorted into two groups that received different doses of rabies virus IM in bilateral masseter muscles. One group received 10^5.9 MICLD₅₀, while the other received 10^6.9 MICLD₅₀ in total volume of 1.0 mL. Rabies virus infection was definitively diagnosed by direct fluorescent antibody assay (DFA) at the Colorado State University Diagnostic Laboratory from collected brainstem and cerebellar tissues [30].

Skunks were individually housed and the day before fecal collection all feces was removed from the pens.Fecal samples that appeared within the last 24 hours were collected from the pen floor, cage catch pan, or den box and the freshest of those fecal samples were preferentially collected [12]. During the vaccination trial, fecal samples were collected at days zero, 29, 60, and 335 post-vaccine (pv). Fecal samples were also collected after the rabies virus challenge on day 14 post-infection, unless the animal was euthanized due to illness before the full 14 days. Three of five animals in the Rinf group, one of three animals given SHvax group, and none of the Rvax animals were euthanized before 14 days post infection. Approximately 2-5g of feces were directly deposited into whirl-pak specimen bag. Samples were stored on ice packs for up to 6 hrs after collection and then stored in -80°C.

Sample processing and sequencing

Samples were extracted using Dneasy® PowerSoil® Kit (Qiagen, Hilden, Germany) according to manufacture directions. We used this extraction kit on all samples to prevent introducing biases from cross-kit comparisons [31]. DNA extracts were stored in -80°C until they could be preppared for sequencing. For each sequencing run we included a commercially available mock community of bacterial DNA, processed through the same library preparation as the samples, for monitoring error and quality control (ZymoBIOMICS microbial community DNA standard, Zymo Research, CA, USA cat. D6306).

To evaluate the fecal microbiome we used a two-step amplicon sequencing approach similar to Galan, Pons [32]. We first amplified the hypervariable region 4 (V4) of the 16S rRNA gene using a multiplex of primers that were based on 515F and 806R primers [33, 34]. We modified the primers to contain different heterogeneity spacers on the 5’ end between partial sequencing adaptors and the amplification primer sequences (S1 Table). All PCR steps were amplified using AccuPrime™ Pfx DNA Polymerase (Invitrogen, Carlsbad, California, USA). The PCR mix contained 25μL of AccuPrime™, 1μL at 10μM for each primer mix, and 5μL of template DNA. The cycling conditions were an initial melting phase at 95°C for 5 min, then 25 cycles of 95°C for 30 sec, 50°C for 60 sec, and 68°C for 98 sec, and a final extension at 68°C for 5 min. To eliminate contamination by primer dimers from 16S amplification, which can decrease sequencing quality, samples were purified using Mag-Bind® TotalPure NGS (Omega Bio-tek, Inc., Norcross, Georgia, USA) beads using a 0.8 volume of beads and eluted in 50μL to dilute samples to improve the second PCR amplification.

The second PCR added sequencing adaptors and dual indices to the purified amplicons (S1 Table). This reaction included a PCR mix containing 25μL of AccuPrime™, 1μL at 10μM of each primer, and 3μL of purified template. The cycling conditions were as follows, 95°C for 5 min, then 8 cycles of 95°C for 30 sec, 50°C for 60 sec, and 68°C for 98 sec, and a final extension at 68°C for 5 min.

As the final purification step, the samples were normalized on SequalPrep™ Normalization Plate Kit (Applied Biosystems, Foster City, California, USA) according to manufacturer’s directions. A subset of samples was checked for normalization of concentration using Qubit™ dsDNA HS Assay Kit (Invitrogen, Carlsbad, California, USA). Samples were then pooled in equal volumes and preppared for paired-end sequencing on an Illumina MiSeq System; we used the 500-cycle MiSeq Reagent Kit v2 (Illumina, San Diego, California, USA).

Bioinformatics

Illumina sequencing fastq files were demultiplexed and read quality was evaluated using FastQC v.0.11.5 [35]. Reads were then quality-filtered using trimmomatic v.36 set at a sliding window of four, optimum PHRED quality score of 25, and a minimum read length of 150 [36]. Files containing zero reads or labeled as unidentified were removed from further data processing to prevent issues in downstream analyses. Taxonomic classification of samples was performed using mothur v.1.44.0 [37] implementing the PE-de-novo-processing script available at the Abdo-Lab/Microbiome-Analysis-Scripts GitHub [37]. Chimera identification within mothur was completed using vsearch v.2.13.3 and silva v.132 as the reference database [38–40].

Statistical analyses

Statistical analyses were completed in R v4.1.3 and RStudio v2022.2.1.461 [41] which included removal of OTUs that were present in extraction blanks (EB) and negative controls. The following statistical analyses were conducted on grouping the sampling timepoints into three groups, Rvax = vaccine treatment, SHvax = sham vaccine, and Rinf = rabies virus infection only. The analyses were conducted on a randomized, repeated measures study design within treatment groups that were self-controlled. Each subgroup included samples collected before any treatment to determine the baseline skunk microbiome within groups (Rvax-start, SHvax-start, Rinf-start), two treatment subgroups (Rvax, sham vaccine SHvax) and three post rabies virus infection subgroups (Rvax-infect, SHvax-infect, Rinf-infect; Table 1, S2 Table). Normality of sequencing depth, number of datapoints, number of OTUs, and number genera per animal was evaluated with a Shapiro-wilks test. We calculated alpha diversity indices, Shannon diversity, Chao1, inverse Simpson (InvSimpson), and rarified species richness, using vegan and phyloseq package v.1.38.0 [42]. We used a rarefaction curve produced by the vegan package v2.6–2 [43] to evaluate the relationship between host-species richness and read depth per sample. Alpha diversity index confirmation to normality was tested using the Shapiro-Wilks test. For indices that were normally distributed, we tested for significant differences among treatment subgroups using a generalized linear model (GLM) and analysis of variance (ANOVA). For indices that deviated from normality we evaluated significance among treatment subgroups using the Kruskal-Wallis rank sum test. Barplots of normalized relative abundance of bacterial genera per treatment subgroup were created using ggplot2 [44] with a minimum relative abundance of 0.01. Venn diagrams of shared OTUs among treatment subgroups were generated with ggvenn v0.1.9 [45] with each group, Rvax, SHvax, and Rinf, shown separately. Visualization of data trends and clustering of microbial community structure per sampling time point (see experimental design in Table 1) was conducted using nonmetric Multidimensional Scaling (NMDS) [46] on OTUs and using Bray-Curtis dissimilarity [46] with OTU data normalized using Cumulative Sum Scaling (CSS) [47]. Significance of beta diversity was assessed among subgroups within each treatment group separately, Rvax, SHvax, and Rinf, using PERMANOVA [48] with 1000 permutations. For significant PERMANOVA we tested for the influence of lack of homogeneity of variance using a beta dispersion test with 1000 permutations. Relative abundance of CSS normalized OTUs was estimated using the phylosmith v1.0.6 package [49]. The Spearman rank correlation coefficient was used to evaluate correlation between OTU relative abundance in the Rvax (vaccination time points only) and fecal VOC peak data from Kimball, Volker [12] to determine if VOC concentration is correlated to specific bacterial taxa in vaccinated animals.

[Figure omitted. See PDF.]

Also presented are the mean (range) of alpha diversity metrics estimated for each treatment subgroup.

We used Random Forests [RF; 50] to determine if any bacterial taxa were predictive of classifying a sample to vaccination state using the Rvax-start, Rvax, and Rvax-infect as subgroups. We first identified the optimal number of features (OTUs) to include in constructing regression trees to be that with the median out-of-bag (OOB) error rate at the asymptote. We used the parameter ntree = 1000 for the total number of trees to grow in the forest and set the parameter importance to TRUE to assess importance of the different features in prediction and ran five iterations to evaluate consistency of assignment error.

Results

Depth, OTUs, Genera, and Alpha diversity

Processing, quality filtering, and removal of samples with fewer than 1000 reads, based on rarefaction curves (S1 Fig), resulted in a total of 53 data points obtained through repeated measurements of 15 animals (seven vaccine (Rvax), three sham vaccine (SHvax), and 5 rabies infection only (Rinf; Table 1)). This resulted in 2,818,526 sequence reads with a mean depth per sample of 53,180 reads (range: 1013–168,770; S2 Fig) and 213 operational taxonomic units (OTUs) representing 75 bacterial genera. Raw sequencing read data are available at the National Center for Biotechnology Information’s (NCBI) Sequence Read Archive (SRA) under accession number PRJNA875681. The relative abundances of each bacterial genera are presented Fig 1. The results of the tests for normality tests for sequencing depth, number of datapoints, number of OTUS, and number of genera are presented in S3 Table, with mean depth and number of datapoints showing significant departures. Sequencing depth, number of OTUs and genera, and values of alpha diversity indices for each treatment group are presented in Table 1 and with mean and median plotted in Fig 2 (depth and alpha diversity for each datapoint is provided in S2 Table with OTUs, and genera in S4 and S5 Tables). Only Shannon diversity conformed to a normal distribution based in the Shapiro-Wilks test (W = 0.98, p = 0.56); Chao1 (W = 0.92, p = 0.0013), richness (W = 0.93, p = 0.0053), and inverse Simpson (W = 0.86, p = 2.2e-05) were not normally distributed. No discernible pattern to the alpha diversity among treatment subgroups is evident from Fig 2 and none of the statistical tests for differences among mean diversities were significant: Chao1 (Kruskal-wallis; χ² = 13.31, p = 0.065), richness (Kruskal-wallis; χ² = 12.06, p = 0.099), inverse Simpson (Kruskal-wallis; χ² = 8.07, p = 0.33), and Shannon (ANOVA; F = 1.14, p = 0.36). We also tested for differences in alpha diversity based on time of sampling and sex of animal and none of the tests were significant. There were five bacterial phlya present across all samples and the dominant phyla in the skunk microbiome were Firmicutes (61% of OTUs) followed by Proteobacteria (15%). The Firmicutes OTUs were represented predominantly by families Clostridiaceae (26%) followed by Lachnospiraceae (16%), Peptostreptococcaceae (14%), Fusobacteriaceae (12%), with other OTUs comprising less than 10% of the total (S2 Table).

[Figure omitted. See PDF.]

Each color represents the proportion of DNA sequences associated with each genus. Each individual bar plot represents a sample datapoint in each of the treatment subgroups. See Table 1, S1 Table, and methods for description of abbreviated treatment subgroups listed on top of plots.

[Figure omitted. See PDF.]

Each plot displays the median as horizontal bar and the mean as a black dot. The metric displayed in each plot is listed on the y-axis and the treatment subgroups is on the x-axis. See Table 1, S1 Table, and methods for definition of treatment subgroups.

Beta diversity

The NMDS plots demonstrate a shift in microbiome diversity among the subgroups within each treatment, Rvax, SHvax, Rinf but each subgroup still showed overlap (Fig 3). However, any structure must be interpreted within the context that only the Rinf plot had stress level below 0.1. But review of the plots demonstrates a shift in microbial communities following rabies virus infection, except in the Rvax group (Fig 3). The PERMANOVA on Bray-Curtis dissimilarity between treatment subgroups was significant for the Rvax group (F = 2.11, p = 0.009) but not the SHvax (F = 0.95, p = 0.49) or Rinf group (F = 1.56, p = 0.11). The shared OTUs among subgroups within each treatment are shown in the Venn diagrams in Fig 4. Notably, the Rvax-infect subgroup only has one unique taxon compared to Rvax-start and Rvax whereas the SHvaxgroup had approximately equal percentages of unique taxa in each subgroup and the Rinf group had a reduction in unique taxa when comparing the Rinf-start and Rinf-infect subgroups. The Rvax group differed from the SHvax group in that the percentage of unique OTUs in the Rvax group was much larger but the taxonomic uniqueness of the Rvax-infect was decreased compared to the SHvax-infect. This unbalance could be due seven sampled animals in the Rvax group versus three animals in the SHvax group.

[Figure omitted. See PDF.]

Individual plots were estimated from Bray-Curtis distance between treatment subgroups. Each subgroup ID represents the group centroid of the 95% confidence ellipsoid.

[Figure omitted. See PDF.]

Each individual plot represents a treatment group. See Table 1, S2 Table, and methods for definition of subgroup abbreviations.

Random forests

Using RF we attempted to identify bacterial OTUs that were important for classifying the samples in the Rvax group into their respective subgroups, Rvax-start, Rvax, and Rvax-infect. Through RF tuning, the median OOB error asymptote was at 44 OTUs with a value of 0.065 (S3 Fig). The mean OOB error rate across 5 runs using 44 OTUs and 1000 trees was 7.1%, which means that approximately 93% of the time samples were assigned to their correct subgroup. The confusion matrix classification error rate for assigning sample to Rvax and Rvax-infect was zero but the mean error rate across the five runs for Rvax-start was 0.34 (S4 Fig). Four to five of the seven Rvax-start samples consistently assigned to the appropriate subgroup. Across four of the RF replicates, two to three of the Rvax-start assigned to Rvax-infect. In a single run, two assigned to Rvax-start and one assigned to Rvax. Fig 5 and S6 Table show the 13 OTUs with the largest influence on cluster assignment. Of these 13, two belong to the family Clostridiaceae, three belong to the Streptococcaceae (specifically the genus Streptococcus), four belong to the family Enterobacteraciae, one belongs to the Erysipelotrichaceae (specifically the genus Turcibacter), two in the Peptostreptococcaceae family, with a single OTU each in the Helicobacteraceae (genus Helicobacter). The important OTUs belong to four orders in three phyla, Proteobacteria, Firmicutes, and Epsilonbacteraeota.

[Figure omitted. See PDF.]

The Gini importance measures on the x-axis represent important OTUs for classifying samples in the vaccine (Rvax) treatment group into subgroups. The y-axis is the OTU and associated taxonomy.

Comparison to VOC data

Spearman rank correlation results for 13 OTUs with the largest mean Gini importance compared to seven VOCs that had a significant positive correlation with vaccination status in skunks from Kimball, Volker [12] are shown in Fig 6 and Table 2. The largest positive correlation was 0.70 and the lowest was -0.57. There were two correlations greater than 0.50 and one less than -0.50. There were 23 instances of p-values less than 0.1 and 14 instances of values less than 0.05 (without correction for multiple comparisons). After correction for multiple comparisons using false discovery rate there were only three significant tests less than 0.1 and all of these were less than 0.05. Two of the significant tests, which were associated with positive correlations, were from comparisons with the Romboustia genus (ninth highest Gini importance of 0.44) and two different VOCs, 2-pentylfuran and unknown compound A. Romboustia also had three more tests that were significant with other VOCs but these were without correction for multiple testing. The remaining significant test was for associations with a Streptococcus OTU and it was a negative correlation coefficient. The unknown compound A was the VOC that had the most correlations with p-values less than 0.1 without correction, and one that was significant following correction.

[Figure omitted. See PDF.]

The size increases and color darkens with increasing or decreasing Spearman rank correlation coefficient. The codes for the VOCs are in S5 Table and values for correlation coefficients are in Table 2. Each of the OTUs was identified in this study as important for classifying samples into vaccination treatment (Rvax) subgroups (see Table 1). The VOCs were identified by Kimball et al., (2019) as correlated with striped skunk vaccination status.

[Figure omitted. See PDF.]

Each of the OTUs was identified in this study as important for classifying samples into vaccination treatment (Rvax) subgroups (see Table 1). The VOCs were identified by Kimball et al,. (2019) as correlated with striped skunk vaccination status and the letters in parentheses are the VOC codes from Kimball et al., 2019. Values in bold were statistically significant (p < 0.10) without correction, * are values that were signficant followig false discovery rate correction, and correlation coefficients over 0.5 and under -0.5 are underlined.

Discussion

The goal of this study was to evaluate whether oral rabies vaccination alters the skunk fecal microbiome and if these changes can be used to predict an animal’s vaccination status. Using samples collected from individuals that were given an oral vaccine, a sham vaccine, and/or infected with rabies virus, we evaluated the alpha and beta diversities, used Random Forests to identify taxa that are indicative of vaccination status, and compared these data to fecal VOC data from Kimball, Volker [12] collected from the same samples. We determined that alpha diversity was not different among the subgroups within the vaccination, sham vaccine, or rabies virus infection treatments. The NMDS plot demonstrated that no major shifts among groups were visible but there were trends. However, the beta diversity was significant between subgroups within the vaccine group (Rvax), but not between subgroups in the sham vaccine or rabies virus infection only groups. Furthermore, there were bacterial taxa that had significant impact on classifying samples into their respective treatment subgroups within the rabies vaccination group. Finally, there were positive correlations between OTUs and VOCs that demonstrated correlation with vaccination status.

Our hypothesis was rooted in previous work that demonstrated microbiomes respond to both vaccination and pathogen infection [e.g., 7, 8, 51]. Counter to these previous studies, alpha diversity did not differ among treatment subgroups. Yet, the beta diversity among the vaccine treatment subgroups was statistically significant. Given that beta diversity was significantly different, and RF identified bacterial taxa that classified of samples into vaccination subgroups, suggest that community composition shifted in response to vaccination. If beta diversity is the only metric that reveals vaccine ingestion, then before and after samples would be required which could limit this method. However, to test the replicability of these results and mechanisms for differences in beta diversity, but not alpha diversity, a more expansive study is needed.

In this study, we documented biodiversity of striped skunk fecal microbiome from a rabies vaccine efficacy trial. The two dominant bacterial phyla identified in our samples were Firmicutes and Proteobacteria. These two phyla have also been identified as prominent in most of the mammalian gastrointestinal tract microbiome studies to date, and specifically, Firmicutes is the dominant phylum in carnivore gut microbiomes [52]. Mammalian gut microbiota are known to be labile and shift in response to changes in diet [53]. For this study, each individual animal was captive raised and fed an identical diet during the experiment and we expect that captive skunk microbiomes may have differences compared to wild skunks, as observed for other mammals [54]. Captive mammals with exposure to humans often demonstrate a “humanization” of the gut microbiome which would impact generalizing inferences from studies on captive animals to wild ones [55]. One aspect that validates the use of captive animals for basic mechanistic responses of microbiome to immunological or physiological changes is that host phylogeny dictates community composition, but not OTU abundance, and mammals may demonstrate the strongest association of phylogeny and microbiota among vertebrate clades [2]. Thus, we can hypothesize that captive skunks would have a core microbiome with broad patterns of community composition and responses similar to wild skunks, however this must be empirically tested before extrapolation. A complete understanding of the effects of phylogeny, life stage, habitat, season, physiological state (e.g., gestation, torpor), and immunological status is critical for validation of a reliable microbiome biomarker in wild animals.

Some, but not all, of the bacteria taxa with the greatest influence on sample classification were correlated with VOCs and/or have been associated with oral vaccination. Streptococcus was the most influential OTU on RF classification which corresponds to previously identified association of gut Streptococcus to inflammation in domestic dogs and cats [56, 57]. Of the highly influential OTUs on RF classification the Enterobacteriaceae were represented among the largest Gini importance factors but were negatively, or slightly positively, correlated with the VOCs hexanal (F), 2-pentylfuran (M), 4-nonanone (R), and unknown compound A (V) of which 2-pentylfuran is a known bacterial metabolite [58]. The Enterobacteriaceae include a wide range of species associated with gut microbiota in mammals with diverse ecological niches. Since these OTUs belong to unclassified genus/genera and are unknown species, little can be said about the role these may play in the response to vaccination and/or infection. Both Abdo, LeCureux [7] and Zhang, Wu [6] found that the order Clostridiales (families Clostridiaceae, Lachnospiraceae (genus Tyzzerella) and Peptostreptococcaceae (genus Romboustia from current study) were associated with vaccination classification and immune response. Specifically, Zhang, Wu [6] identified that abundance of unclassified Clostridiales and Lachnospiraceae OTUs were positively correlated with rabies vaccine immune response in mice.

We found Romboustia had the strongest correlation with VOCs and other studies corroborate our results that this genus is associated with metabolites dependent on disease and study model, and with dysbiosis in the microbiome and host inflammation. In a murine ApoE^-/^- model for atherosclerosis, Yan, Wang [59] determined metformin treatment altered microbiome diversity in a treatment-dependent manner; specifically, Romboutsia decreased in abundance due to metformin treatment and was positively associated with inflammatory immune markers and negatively correlated with various short chain fatty acids. Jin, Jia [60] also found Romboutsia was positively correlated with TNFα and the metabolites methoxy-4-hydroxypheylglycol sulfate and cis-(6,9,12)-Linoleic acid in Bamei suckling piglets given a milk replacer supplement, and negatively correlated with genes regulating intestinal barrier permeability and function (including ZO-1) however association with specific immune functions has not been determined.

Logical next steps to improve and expand this study would be to increase treatment group sample sizes, evaluate wild-caught striped skunks versus captive raised, test other rabies reservoir species and targets of management such as raccoons, and determine if the ingestion of the oral bait to deliver the rabies vaccine versus direct instillation has similar effects on the host microbiome response. The oral rabies vaccine bait used for raccoon and canid rabies management in the U.S. uses tetracycline as a biomarker. Detecting the biomarker is invasive as it requires anesthetizing an animal and removing a tooth [29]. Identifying a fecal biomarker of vaccine bait uptake would improve the ease of sample collection for population monitoring and reduce the risk to field biologists and animal welfare, but characterization of the limits to sample freshness, number of samples, and regularity of collection would all need to be investigated. Another interesting test would be to explore the impact of tetracycline in the oral bait on the gut microbiome. Tetracycline is an antibiotic and would likely impact the gut microbiota [61]. Perhaps in carnivores the ingestion of tetracycline can leave a signature in the fecal microbiome that would allow for detection or recent bait ingestion. The oral vaccine used in this study did not contain tetracycline but the MEM culture media used to dilute the vaccine and prepare sham treatments did contain low levels of an antibiotic and antimycotic. But our study design accounted for this by incorporating the sham vaccine to test for artifacts of either animal handling or the vaccine media. We saw different responses in the vaccinated versus sham groups which suggests that the MEM media was not a major factor the experimental outcomes.

Biomarkers have high potential to facilitate non-invasive monitoring of populations for pathogens and vaccination status, which is critical to wildlife disease management. Continually pushing the boundaries of sample types and diagnostic molecules will lead to increased efficiency, decreased cost and risk, and open new doors to additional types of diagnostic samples for wildlife disease surveillance. In this study we identified some correlations between fecal bacterial taxa, vaccination status, and VOCs. While we cannot infer beyond the results to identify causal mechanisms, these results can lead to additional hypotheses and further studies about new biomarkers to assist with rabies, and other pathogen surveillance in wildlife.

Supporting information

S1 Table. Indexes.

Indexes used for each sample in this study of skunk microbiome.

https://doi.org/10.1371/journal.pone.0285852.s001

(XLSX)

S2 Table. Sample metadata and diversity metrics.

Summary data for each sample datapoint (ID) collected for the skunk microbiome study. Includes are Animal ID, animal sex, date of collection, dose of rabies virus infection, the treatment group associates with data analysis (Table 1), treatment subgroup, days post vaccination, days post infection, rabies virus neutralizing antibody (rVNA) from Kimball et al., 2019, and alpha diversity metrics based on OTUS for each individual sampling datapoint.

https://doi.org/10.1371/journal.pone.0285852.s002

(XLSX)

S3 Table. Shapiro-Wilks test.

Results from Shapiro-Wilks test for normality of the data across all samples.

https://doi.org/10.1371/journal.pone.0285852.s003

(XLSX)

S4 Table. OTUs and taxonomy.

Operational taxonomic units (OTU) and associated taxonomy of 16S bacterial sequences collected from striped skunk fecal samples.

https://doi.org/10.1371/journal.pone.0285852.s004

(XLSX)

S5 Table. Sequencing depth.

Operational taxonomic units (OTUs) of bacterial 16S DNA sequences collected from striped skunk fecal samples and sequence depth for each sample data point. See S1 Table for treatment groups associated with datapoints.

https://doi.org/10.1371/journal.pone.0285852.s005

(XLSX)

S6 Table. Random Forests OUT importance.

Operational taxonomic units (OTUs) and associated taxonomy identified from bacterial 16S DNA sequences collected from striped skunk fecal samples. The relative mean Gini importance (a) estimated by random forests for 44 OTUs and 1000 trees that influence the classification of TV datapoints into treatment subgroups: TV-start, TV, and TV-infect.

https://doi.org/10.1371/journal.pone.0285852.s006

(XLSX)

S7 Table. VOCs from Kimball et al. (2019).

Volatile organic compounds (VOC) identified in striped skunk fecal samples and their code used in Fig 6. Table was adapted from Kimball et al., (2019).

https://doi.org/10.1371/journal.pone.0285852.s007

(XLSX)

S1 Fig. Rarefaction curves based on operational taxonomic units from skunk fecal microbiome.

https://doi.org/10.1371/journal.pone.0285852.s008

S2 Fig. Histogram of sequencing depth per sample for skunk fecal microbiome samples.

https://doi.org/10.1371/journal.pone.0285852.s009

S3 Fig. Random Forests Tuning results for skunk fecal microbometaxa based on 1000 trees.

X-axis Represents the number of OTUS included in the permutation and y-axis is the median out of bag error.

https://doi.org/10.1371/journal.pone.0285852.s010

S4 Fig. Random Forests classification error rate for skunk fecal microbiome samples representing vaccinated (TV), Before any treatment (TV-start), and post-vaccinrabies virus challenge (TV-infect).

The lines represent the out of bag (OOB) error rate (y-axis) for the number of trees included in the analysis (x-axis).

https://doi.org/10.1371/journal.pone.0285852.s011

Acknowledgments

The authors thank the following persons for invaluable technical assistance with the prior live animal work: N. Walker, S. Eaton, J. Kanine, C. Wickham, D. Wostenberg, D. Kohler, C. Ellis, S. Volker, and the Animal Care Unit at NWRC. We thank A. Davis for guidance through the writing of this manuscript. The authors would like to thank and recognize the USDA APHIS Science Fellows Program for making this work possible.

Citation: Hopken MW, Gilfillan D, Gilbert AT, Piaggio AJ, Hilton MS, Pierce J, et al. (2023) Biodiversity indices and Random Forests reveal the potential for striped skunk (Mephitis mephitis) fecal microbial communities to function as a biomarker for oral rabies vaccination. PLoS ONE 18(8): e0285852. https://doi.org/10.1371/journal.pone.0285852

About the Authors:

Matthew W. Hopken

Roles: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

E-mail: [email protected]

Affiliations: United States Department of Agriculture, Animal and Plant Health Inspection Service, Wildlife Services, National Wildlife Research Center, Fort Collins, Colorado, United States of America, Department of Microbiology, Immunology, and Pathology, Colorado State University, Fort Collins, Colorado, United States of America

ORICD: https://orcid.org/0000-0003-3861-6153

Darby Gilfillan

Roles: Data curation, Formal analysis, Methodology, Writing – original draft, Writing – review & editing

Affiliation: Department of Microbiology, Immunology, and Pathology, Colorado State University, Fort Collins, Colorado, United States of America

Amy T. Gilbert

Roles: Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – review & editing

Affiliation: United States Department of Agriculture, Animal and Plant Health Inspection Service, Wildlife Services, National Wildlife Research Center, Fort Collins, Colorado, United States of America

Antoinette J. Piaggio

Roles: Conceptualization, Funding acquisition, Resources, Writing – review & editing

Affiliation: United States Department of Agriculture, Animal and Plant Health Inspection Service, Wildlife Services, National Wildlife Research Center, Fort Collins, Colorado, United States of America

Mikaela Samsel Hilton

Roles: Investigation, Methodology, Writing – review & editing

Affiliation: United States Department of Agriculture, Animal and Plant Health Inspection Service, Wildlife Services, National Wildlife Research Center, Fort Collins, Colorado, United States of America

James Pierce

Roles: Investigation, Methodology, Writing – review & editing

Affiliation: Department of Microbiology, Immunology, and Pathology, Colorado State University, Fort Collins, Colorado, United States of America

Bruce Kimball

Roles: Conceptualization, Resources, Writing – review & editing

Affiliation: Monell Chemical Senses Center, Philadelphia, Pennsylvania, United States of America

Zaid Abdo

Roles: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

Affiliation: Department of Microbiology, Immunology, and Pathology, Colorado State University, Fort Collins, Colorado, United States of America

References

1. Hird SM. Evolutionary Biology Needs Wild Microbiomes. Frontiers in Microbiology. 2017;8. pmid:28487687

2. Youngblut ND, Reischer GH, Walters W, Schuster N, Walzer C, Stalder G, et al. Host diet and evolutionary history explain different aspects of gut microbiome diversity among vertebrate clades. Nature Communications. 2019;10(1):2200. pmid:31097702

3. Williams CL, Caraballo-Rodríguez AM, Allaband C, Zarrinpar A, Knight R, Gauglitz JM. Wildlife-microbiome interactions and disease: exploring opportunities for disease mitigation across ecological scales. Drug Discovery Today: Disease Models. 2018;28:105–15. https://doi.org/10.1016/j.ddmod.2019.08.012.

4. Trevelline BK, Fontaine SS, Hartup BK, Kohl KD. Conservation biology needs a microbial renaissance: a call for the consideration of host-associated microbiota in wildlife management practices. Proceedings of the Royal Society B: Biological Sciences. 2019;286(1895):20182448. pmid:30963956

5. Desselberger U. The Mammalian Intestinal Microbiome: Composition, Interaction with the Immune System, Significance for Vaccine Efficacy, and Potential for Disease Therapy. Pathogens. 2018;7(3):57. pmid:29933546

6. Zhang Y, Wu Q, Zhou M, Luo Z, Lv L, Pei J, et al. Composition of the murine gut microbiome impacts humoral immunity induced by rabies vaccines. Clinical and Translational Medicine. 2020;10(4):e161. pmid:32898335

7. Abdo Z, LeCureux J, LaVoy A, Eklund B, Ryan EP, Dean GA. Impact of oral probiotic Lactobacillus acidophilus vaccine strains on the immune response and gut microbiome of mice. PLOS ONE. 2019;14(12):e0225842. pmid:31830087

8. Guevarra RB, Cho JH, Cho JH, Lee JH, Kim H, Kim S, et al. Oral Vaccination against Lawsoniaintracellularis Changes the Intestinal Microbiome in Weaned Piglets. Animals. 2021;11(7):2082. pmid:34359211

9. Lynn MA, Tumes DJ, Choo JM, Sribnaia A, Blake SJ, Leong LEX, et al. Early-Life Antibiotic-Driven Dysbiosis Leads to Dysregulated Vaccine Immune Responses in Mice. Cell Host & Microbe. 2018;23(5):653–60.e5. pmid:29746836

10. Ellis CK, Stahl RS, Nol P, Waters WR, Palmer MV, Rhyan JC, et al. A Pilot Study Exploring the Use of Breath Analysis to Differentiate Healthy Cattle from Cattle Experimentally Infected with Mycobacterium bovis. PLOS ONE. 2014;9(2):e89280. pmid:24586655

11. Stahl RS, Ellis CK, Nol P, Waters WR, Palmer M, VerCauteren KC. Fecal Volatile Organic Ccompound Profiles from White-Tailed Deer (Odocoileus virginianus) as Indicators of Mycobacterium bovis Exposure or Mycobacterium bovis Bacille Calmette-Guerin (BCG) Vaccination. PLOS ONE. 2015;10(6):e0129740. pmid:26060998

12. Kimball BA, Volker SF, Griffin DL, Johnson SR, Gilbert AT. Volatile metabolomic signatures of rabies immunization in two mesocarnivore species. PLOS Neglected Tropical Diseases. 2019;13(12):e0007911. pmid:31790413

13. Kimball BA, Yamazaki K, Kohler D, Bowen RA, Muth JP, Opiekun M, et al. Avian Influenza Infection Alters Fecal Odor in Mallards. PLOS ONE. 2013;8(10):e75411. pmid:24146753

14. Gomez A, Petrzelkova K, Yeoman CJ, Vlckova K, Mrázek J, Koppova I, et al. Gut microbiome composition and metabolomic profiles of wild western lowland gorillas (Gorilla gorilla gorilla) reflect host ecology. Molecular Ecology. 2015;24(10):2551–65. pmid:25846719

15. Carey HV, Walters WA, Knight R. Seasonal restructuring of the ground squirrel gut microbiota over the annual hibernation cycle. American Journal of Physiology-Regulatory, Integrative and Comparative Physiology. 2013;304(1):R33–R42. pmid:23152108.

16. Shi Z, Zou J, Zhang Z, Zhao X, Noriega J, Zhang B, et al. Segmented Filamentous Bacteria Prevent and Cure Rotavirus Infection. Cell. 2019;179(3):644–58.e13. pmid:31607511

17. Ganz HH, Doroud L, Firl AJ, Hird SM, Eisen JA, Boyce WM, et al. Community-Level Differences in the Microbiome of Healthy Wild Mallards and Those Infected by Influenza A Viruses. mSystems. 2017;2(1):e00188–16. pmid:28293681

18. Taberlet P, Luikart G. Non-invasive genetic sampling and individual identification. Biological Journal of the Linnean Society. 2008;68(1–2):41–55.

19. Nel LH, Markotter W. Lyssaviruses. Critical Reviews in Microbiology. 2007;33(4):301–24. pmid:18033596

20. Troupin C, Dacheux L, Tanguy M, Sabeta C, Blanc H, Bouchier C, et al. Large-Scale Phylogenomic Analysis Reveals the Complex Evolutionary History of Rabies Virus in Multiple Carnivore Hosts. PLOS Pathogens. 2016;12(12):e1006041. pmid:27977811

21. Kuzmin IV, Shi M, Orciari LA, Yager PA, Velasco-Villa A, Kuzmina NA, et al. Molecular Inferences Suggest Multiple Host Shifts of Rabies Viruses from Bats to Mesocarnivores in Arizona during 2001–2009. PLOS Pathogens. 2012;8(6):e1002786. pmid:22737076

22. Rupprecht CE, Glickman LT, Spencer PA, Wiktor TJ. EPIDEMIOLOGY OF RABIES VIRUS VARIANTS: DIFFERENTIATION USD MONOCLONAL ANTIBODIES AND DISCRIMINANT ANALYSIS. American Journal of Epidemiology. 1987;126(2):298–309. pmid:3300280

23. Ma X, Bonaparte S, Toro M, Orciari LA, Gigante CM, Kirby JD, et al. Rabies surveillance in the United States during 2020. Journal of the American Veterinary Medical Association. 2022;260(10):1157–65. pmid:35522584

24. Slate D, Algeo TP, Nelson KM, Chipman RB, Donovan D, Blanton JD, et al. Oral Rabies Vaccination in North America: Opportunities, Complexities, and Challenges. PLOS Neglected Tropical Diseases. 2009;3(12):e549. pmid:20027214

25. Elmore SA, Chipman RB, Slate D, Huyvaert KP, VerCauteren KC, Gilbert AT. Management and modeling approaches for controlling raccoon rabies: The road to elimination. PLOS Neglected Tropical Diseases. 2017;11(3):e0005249. pmid:28301480

26. Guerra MA, Curns AT, Rupprecht CE, Hanlon CA, Krebs JW, Childs JE. Skunk and raccoon rabies in the eastern United States: temporal and spatial analysis. Emerging Infectious Diseases. 2003;9(9):1143. pmid:14519253

27. Wallace RM, Gilbert A, Slate D, Chipman R, Singh A, Cassie W, et al. Right Place, Wrong Species: A 20-Year Review of Rabies Virus Cross Species Transmission among Terrestrial Mammals in the United States. PLOS ONE. 2014;9(10):e107539. pmid:25295750

28. Davis AJ, Kirby JD, Chipman RB, Nelson KM, Gilbert AT. Data-Driven Management—A Dynamic Occupancy Approach to Enhanced Rabies Surveillance Prioritization. Viruses. 2021;13(9):1795. pmid:34578376

29. Algeo TP, Norhenberg G, Hale R, Montoney A, Chipman RB, Slate D. ORAL RABIES VACCINATION VARIATION IN TETRACYCLINE BIOMARKING AMONG OHIO RACCOONS. Journal of Wildlife Diseases. 2013;49(2):332–7. pmid:23568908

30. Gilbert A, Johnson S, Walker N, Beath A, VerCauteren K. FLAVOR PREFERENCE AND EFFICACY OF VARIABLE DOSE ONTARIO RABIES VACCINE BAIT (ONRAB) DELIVERY IN STRIPED SKUNKS (MEPHITIS MEPHITIS). Journal of Wildlife Diseases. 2018;54(1):122–32. pmid:29077542

31. Kennedy NA, Walker AW, Berry SH, Duncan SH, Farquarson FM, Louis P, et al. The Impact of Different DNA Extraction Kits and Laboratories upon the Assessment of Human Gut Microbiota Composition by 16S rRNA Gene Sequencing. PLOS ONE. 2014;9(2):e88982. pmid:24586470

32. Galan M, Pons J-B, Tournayre O, Pierre É, Leuchtmann M, Pontier D, et al. Metabarcoding for the parallel identification of several hundred predators and their prey: Application to bat species diet analysis. Molecular Ecology Resources. 2018;18(3):474–89. pmid:29288544

33. Apprill A, McNally S, Parsons R, Weber L. Minor revision to V4 region SSU rRNA 806R gene primer greatly increases detection of SAR11 bacterioplankton. Aquatic Microbial Ecology. 2015;75(2):129–37.

34. Parada AE, Needham DM, Fuhrman JA. Every base matters: assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environmental Microbiology. 2016;18(5):1403–14. pmid:26271760

35. Andrews S. A Quality Control Tool for High Throughput Sequence Data 2010. Available from: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

36. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. pmid:24695404

37. 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. Applied and Environmental Microbiology. 2009;75(23):7537–41. pmid:19801464

38. Rognes T, Flouri T, Nichols B, Quince C, Mahé F. VSEARCH: a versatile open source tool for metagenomics. PeerJ. 2016;4:e2584. pmid:27781170

39. 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 Research. 2007;35(21):7188–96. pmid:17947321

40. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Research. 2012;41(D1):D590–D6. pmid:23193283

41. R_Core_Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria: R Core Team; 2022.

42. McMurdie PJ, Holmes S. phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data. PLOS ONE. 2013;8(4):e61217. pmid:23630581

43. Oksanen J, Simpson GL, Blanchet G, Kindt R, Legendre P, Minchin PR, et al. vegan: Community Ecology Package https://CRAN.R-project.org/package=vegan2022. Available from: https://CRAN.R-project.org/package=vegan.

44. Villanueva RAM, Chen ZJ. ggplot2: elegant graphics for data analysis. Taylor & Francis; 2019.

45. Linlin Y. ggvenn: Draw Venn Diagram by ’ggplot2’ https://CRAN.R-project.org/package=ggvenn2021. Available from: https://CRAN.R-project.org/package=ggvenn.

46. Legendre P, Legendre L. Numerical Ecology. 2nd ed. Amsterdam, The Netherlands: Elsevier B.V.; 1998.

47. Paulson JN, Stine OC, Bravo HC, Pop M. Differential abundance analysis for microbial marker-gene surveys. Nature Methods. 2013;10(12):1200–2. pmid:24076764

48. Anderson MJ. A new method for non-parametric multivariate analysis of variance. Austral Ecology. 2001;26(1):32–46. https://doi.org/10.1111/j.1442-9993.2001.01070.pp.x.

49. Smith SD. phylosmith: an R-package for reproducible and efficient microbiome analysis with phyloseq-objects https://schuyler-smith.github.io/phylosmith/2019. Available from: https://schuyler-smith.github.io/phylosmith/.

50. Breiman L. Random Forests. Machine Learning. 2001;45(1):5–32.

51. Fox BE, Vilander AC, Gilfillan D, Dean GA, Abdo Z. Oral Vaccination Using a Probiotic Vaccine Platform Combined with Prebiotics Impacts Immune Response and the Microbiome. Vaccines. 2022;10(9):1465. pmid:36146543

52. Nishida AH, Ochman H. Rates of gut microbiome divergence in mammals. Molecular Ecology. 2018;27(8):1884–97. pmid:29290090

53. Muegge BD, Kuczynski J, Knights D, Clemente JC, González A, Fontana L, et al. Diet Drives Convergence in Gut Microbiome Functions Across Mammalian Phylogeny and Within Humans. Science. 2011;332(6032):970–4. pmid:21596990

54. McKenzie VJ, Song SJ, Delsuc F, Prest TL, Oliverio AM, Korpita TM, et al. The Effects of Captivity on the Mammalian Gut Microbiome. Integrative and Comparative Biology. 2017;57(4):690–704. pmid:28985326

55. Trevelline BK, Moeller AH. Robustness of Mammalian Gut Microbiota to Humanization in Captivity. Frontiers in Ecology and Evolution. 2022;9. pmid:35291481

56. Pilla R, Suchodolski JS. The Role of the Canine Gut Microbiome and Metabolome in Health and Gastrointestinal Disease. Front Vet Sci. 2020;6:498–. pmid:31993446.

57. Hashimoto-Hill S, Alenghat T. Inflammation-Associated Microbiota Composition Across Domestic Animals. Frontiers in Genetics. 2021;12. pmid:34239536

58. NCBI. PubChem Compound Summary for CID 19602, 2-Pentylfuran. 2022 [September 14, 2022]. Available from: https://pubchem.ncbi.nlm.nih.gov/compound/2-Pentylfuran.

59. Yan N, Wang L, Li Y, Wang T, Yang L, Yan R, et al. Metformin intervention ameliorates AS in ApoE-/- mice through restoring gut dysbiosis and anti-inflammation. PLOS ONE. 2021;16(7):e0254321. pmid:34264978

60. Jin J, Jia J, Zhang L, Chen Q, Zhang X, Sun W, et al. Jejunal inflammatory cytokines, barrier proteins and microbiome-metabolome responses to early supplementary feeding of Bamei suckling piglets. BMC Microbiology. 2020;20(1):169. pmid:32552686

61. Chopra I, Roberts M. Tetracycline Antibiotics: Mode of Action, Applications, Molecular Biology, and Epidemiology of Bacterial Resistance. Microbiology and Molecular Biology Reviews. 2001;65(2):232–60. pmid:11381101