Scat samples were collected during annual trap monitoring of U. littoralis during the fall and winter of 2014–2015 by managers or contractors on the Channel Islands off the coast of southern California, United States (San Miguel Island, SMI; Santa Rosa Island, SRI; Santa Cruz Island, SCZ; Santa Catalina Island, CAT; San Clemente Island, SCL; San Nicolas Island, SNI; Figure 1). See Supplementary Table 1 for sample details. Samples were collected from traps of individual foxes, so the identity of the foxes were known. Scat samples from captive U. littoralis were collected opportunistically in 2014 (Orange County Zoo, OCZ; Santa Barbara Zoo, SBZ). Scat samples were collected in ethanol and then stored at −80°C until DNA extraction could be performed. Total genomic DNA was extracted from the scat samples using the E.Z.N.A Genomic DNA Isolation Kit for stool samples (Omega Bio-Tek, Norcross, GA, United States) following the manufacturer’s instructions with minor changes. Modifications to the protocol included incubating with Proteinase K for 20 min (an increase of 10 min) and eluting the DNA with 100 μL of elution buffer. Glass beads were used during the extractions to ensure bacteria were lysed, per the extraction kit’s instructions. Each batch was extracted alongside an experimental blank containing molecular grade water instead of sample to monitor potential contamination, and samples were extracted no more than 8 months after being collected in the field. Genomic DNA was quantified using an Invitrogen Qubit 2.0 Fluorometer (Carlsbad, CA, United States).

The V4 hypervariable region of the 16S rRNA gene was amplified from DNA extracted from U. littoralis scat using bacterial “universal” primers 515F and 806R, to which 4 ambiguous bases and variable barcodes on the forward primer and variable indexes on the reverse primers were added (Supplementary Table 2; Bartram et al., 2011; Caporaso et al., 2011). We followed a modified Earth Microbiome Project protocol (Gilbert et al., 2014). Briefly, the 16S V4 region was amplified in triplicate reactions using Phusion High-Fidelity DNA Polymerase (New England Biolabs, Ipswich, MA, United States) and then the reactions were cleaned with AMPure XP beads (Beckman Coulter, Brea, CA, United States) at a 1:1 ratio of beads to PCR product volume. The amplicon quantity was again evaluated with the Qubit and the quality was further evaluated with a Bioanalyzer analysis (Agilent Technologies, Santa Clara, CA, United States). Finally, samples were pooled equimolarly for sequencing. We sequenced 95 samples (Supplementary Table 1) as 300 bp paired-end reads on an Illumina MiSeq at UCLA’s Genotyping Core (Los Angeles, CA, United States). Blank samples (N = 5), technical replicates (repeat sequencing library preparations on the same DNA extraction; N = 10), and a positive control (an even mock community) were sequenced as part of the 95 samples. For the Orange County Zoo samples (sample 1 = OCZ.1a and OCZ.1b) and one Santa Barbara Zoo sample (sample 102003 = SBZ.2a and SBZ.2b) biological replicates were sequenced (meaning two different scat samples from the same organism). These two zoo individuals (SBZ.2a/b and OCZ.1a/b) were born in captivity but are descended from San Clemente Island foxes. One zoo sample (SBZ.1) was born on San Clemente Island but was brought into captivity.

Leading ambiguous bases were trimmed from forward reads with Trimmomatic v. 0.33 (Bolger et al., 2014). Barcodes were extracted with extract_barcodes.py, and samples were demultiplexed in QIIME1 v. 1.9.1 (Caporaso et al., 2010). Demultiplexed reads were merged into one artifact with QIIME1 and then imported back into QIIME2. Sequences were processed using denoise-paired in the DADA2 v. 2018.6.0 plugin (Callahan et al., 2016) in QIIME2. Processing included the following steps with default settings unless otherwise noted: denoising, trimming primers (–p-trim-left-f 19, –p-trim-left-r 20), removing chimeras and dereplicate reads, trimming forward reads to a minimum of 270 bp and reverse reads to a minimum of 210 bp based on quality scores, filtering out targets below those length cutoffs, and inferring amplicon sequence variants (ASVs) for both forward and reverse reads. Non-clustering methods that generate sequence variants are convincingly argued to be preferred over traditional 97% sequence clustering approaches due to their higher resolution and reproducibility (Callahan et al., 2017). Taxonomy was assigned using feature-classifier classify-sklearn (Pedregosa et al., 2011) in QIIME2 using the Greengenes database v. 13-8 (DeSantis et al., 2006), then non-target ASVs (mitochondria, archaea, and chloroplasts) were removed. The mock community sample was removed from the dataset for analysis.

The data were imported into RStudio v. 1.1.456 (RStudio Team, 2018) with R v. 3.5.1 (R Core Team, 2016) using the Qiime2r v. 0.99.1 package (Bisanz, 2018) creating an object to use with the Bioconductor R package Phyloseq v. 1.24.2 (McMurdie and Holmes, 2013). Extraction blanks were addressed first by identifying potential contaminants using Decontam v. 1.0.0 (Davis et al., 2018). We used the “combined” method, which is based on frequency (variation based on DNA concentration) and prevalence (increased prevalence in negative controls) probabilities in tandem with a Fisher’s combined probabilities test. This method identified 12 ASVs as potential contaminants which were then removed from the dataset (Supplementary Table 3). After decontamination, extraction blanks were also removed from the dataset. Then technical replicates were merged using mean counts in Phyloseq [see Supplementary Figure 1 for principal coordinates analysis (PCoA) of technical replicates]. The ASV count table was transformed using a variance stabilizing transformation (Anders and Huber, 2010; McMurdie and Holmes, 2014) in DESeq2 v. 1.20.0 (type = “poscounts” to estimate size factors; Love et al., 2014) for beta diversity assessment. For use in Phyloseq (McMurdie and Holmes, 2013), negative transformed values were set to zero as advised by the authors¹.

The dominant taxa present in the scat bacterial microbiome were identified based on ASV proportions (not transformed), and differences among sources were evaluated with analysis of variance (ANOVA) and Tukey HSD post hoc tests. First we agglomerated the count table at the phylum level for all phyla except for the Proteobacteria phylum which was split into its major classes for analysis due to its complexity (Hug et al., 2016), and phyla representing less than 5% of the ASVs were grouped together as “other.” With similar methods, we looked at the top genera with the largest proportions. To address beta diversity, a Bray–Curtis dissimilarity matrix was calculated on variance-stabilizing transformed data and used with hclust (method = ward.D2) to generate a hierarchical clustering dendrogram. The transformed data were also used to conduct a principal coordinates analysis (PCoA) based on Bray–Curtis dissimilarities with Phyloseq.

Testing for statistical differences in scat bacterial communities was done on the transformed ASV counts with Vegan v. 2.5-3 (Oksanen et al., 2020). Permutational ANOVAs were calculated using adonis (permutations = 99999) based on Bray–Curtis dissimilarity (Anderson, 2006, 2017) only on samples with no missing data (Supplementary Table 4A). The data were transformed separately for each dataset after samples with missing data were removed. We did this across all samples and then again for each island and the captive versus wild dataset. The assumption of homogeneity in dispersion within groups was tested with betadisper. We evaluated the best models by including all variables regardless of dispersion and then we sequentially removed variables with low variance explained (R²). The models were evaluated with Akaike information criterion (AIC), taking the lowest value as the best model that was used in the permutational ANOVA. The source code for the AIC method is available at Kdyson (2020). For the combined dataset that included all islands, we tested the effects of age, body condition, month and year collected, sex, source (each island or zoo), weight, and the extraction group in which the sample belonged. For each island, we evaluated the same variables excluding source and year collected when it did not vary from month collected (SNI). For the captive dataset, we included the following variables: captivity status, age, month and year collected, sex, and if the individual was born in captivity or born on an island then taken into captivity. Body condition is a common measurement scale used to gauge animal health (Coonan et al., 2014) that goes from emaciated (body condition = 1) to obese (body condition = 5) and was recorded by wildlife managers in the field. We then examined differential abundances in ASVs for comparisons with significant permutational ANOVA results using the R package DESeq2 and a significance level of 0.01 (with the input being the untransformed, but filtered, count table). We also looked at differentially abundant genera by first agglomerating the count table at the genus level. We looked at alpha diversity in comparisons with significant permutational ANOVA results by calculating Chao1 richness index and Shannon diversity index using untransformed count tables and plot_richness in Phyloseq to understand relative differences.

Bacteroidetes, Firmicutes, Fusobacteria, and Proteobacteria (classes Gammaproteobacteria and Betaproteobacteria) were the dominant phyla in the scat samples; the first two accounting for over 67% of the recovered sequences and the top 6 major taxa collectively accounting for about 94% of the total sequences in each group (Figure 2A). Seven major taxa were significantly different between sources. Firmicutes had the highest number of significant pairwise comparisons and was mainly elevated in zoo foxes compared to island foxes (Supplementary Table 5A and Supplementary Figure 2). Foxes from SMI and SNI were significantly different between 3 out of 7 major taxa, the most of any pairwise comparison (Supplementary Table 5A). Largely foxes from SMI had significantly elevated proportions of Fusobacteria (Supplementary Table 5A and Supplementary Figure 2). Bacteroidia made up the largest average proportion (37.8 ± 12.2%; mean ± 1SD) of identified classes followed by Clostridia (27.4 ± 10.0%). The same pattern held for order; Bacteroidales and Clostridiales made up the largest represented orders at comparable percentages to their respective classes. The top genera were Bacteroides and Fusobacterium accounting for over 38% of the recovered sequences (Supplementary Table 5B and Figure 2B). Of the top genera, 6 were significantly different between sources. Escherichia was elevated in foxes from CAT compared to SRI and the other southern islands of SCL and SNI. Most notably, however, Clostridium and Streptococcus were significantly elevated in captive foxes from OCZ compared to every other source (Supplementary Table 5B and Figure 2B).

The bacterial communities largely clustered into two main groups based on geography (Figures 3A,B). We identified two main clusters based on the dendrogram (Figure 3A). One cluster is entirely composed of samples from the southern islands (SNI, SCL, and CAT), while the other contained all of the northern and zoo samples, as well as a subset of southern samples including about half of the SCL samples and two from SNI. In the PCoA (Figure 3B), axis 1 (16%) largely described the distance between the northern/zoo samples group and the southern islands group, including the aforementioned exceptions, and axis 2 (11.4%) predominantly described the variation within those groups.

The strongest driver of differences among bacterial communities was source (Table 1 and Figure 3C). Across all samples, weight and month of sample collection also significantly influenced the bacterial communities (Table 1). Although, it is important to note that there is a known significant difference in body size among the islands (Wayne et al., 1991), which was also found in our weight data (ANOVA_source, P = 5.0 × 10^–9). The highest mean weight was found on CAT, and SNI had the lowest mean weight. When broken out by island, month was the driving factor on SRI, sex the largest factor on CAT, and on SCL weight and year drove the difference in bacterial communities (Table 1 and Supplementary Figure 3). To look further into the influence of collection timing on SRI and SCL, we evaluated precipitation over the period assessed for each island (Supplementary Figure 4). Sex was also found to be the driving factor on SNI; however, the pattern was driven by one male, and when removed there were no significant factors (Table 1 and Supplementary Figure 3). We investigated the difference between the southern samples that grouped with the northern samples and those that did not and found the best model included weight, month collected, and the group in which the sample was extracted as significant drivers of bacterial community difference. For results from our model testing see Supplementary Tables 4B–H. It should be noted that not all variables were homogenously distributed, but the heterogeneity in the source term is likely biologically relevant (Supplementary Table 6 and Figure 3C). Further, SCZ did not have enough samples (N = 3) to conduct the permutational ANOVA, and there were no significant factors identified for SMI likely due to small sample size (N = 6) after removing samples with missing data.

We found 173 differentially abundant ASVs significant at the α = 0.01 level in our significant factors across islands. On SRI we found 72% of the significantly differentially abundant ASVs were elevated in samples collected in November, while 28% of the ASVs were elevated in December (Table 2). On CAT we found 38% of the significantly differentially abundant ASVs were elevated in female samples, opposed to 62% of the ASVs being elevated in males (Table 2). Additionally, on SCL we found the majority of the significantly differentially abundant ASVs were elevated in samples collected from 2014 (87%) and heavier foxes (88%). Since weight was evaluated as a continuous variable, this result can be interpreted as an increase in ASV abundance with increasing weight. Both on SRI and SCL the genus Megamonas was significantly differentially abundant in samples collected in December and January (for additional differences at the genus level see Supplementary Table 7).

To explicitly test for the effect of captivity on the scat bacterial microbiome, we probed differences between the captive and SCL samples because all zoo individuals were either descendants of SCL foxes or born on the island and then taken into captivity (Figure 4). We found that the best model for these data included the captivity factor along with the month and year that the sample was collected, and all three terms were significant (Table 1). Body condition and weight could not be considered due to missing data. We found that the majority of the significantly differentially abundant ASVs were elevated in wild samples (81%) compared to captive (Table 2). There was a nearly equal number of significantly differentially abundant ASVs between sampled years, but more differentiated ASVs from samples collected in months later in the year. In fact, Megamonas was again seen as a significant player in samples from later months, which could be driven by the same result found in SCL samples alone (Supplementary Table 7).

Statistical analyses on alpha diversity can be misleading (Willis, 2019), therefore we explore qualitative alpha diversity patterns to gain insight into the statistically significant factors we identified through the permutational multivariate analysis of variance (PERMANOVAs). When comparing across sources, scat bacterial communities from the three southern islands, especially CAT, generally had higher ASV richness and evenness compared to the northern Channel Islands (Supplementary Figure 5A). Across all samples analyzed, ASV richness and evenness peaked in October, while richness, but not evenness, trended toward increasing with weight (Supplementary Figures 5B,C). In bacterial communities from SRI, scat collected in November had higher richness and evenness compared to those collected in December (Supplementary Figure 5D). Males had higher evenness, but not richness, compared to females on CAT (Supplementary Figure 5E). Meanwhile, the evenness of bacterial communities on SCL declined with increasing weight while richness was constant (Supplementary Figure 5F). Both richness and evenness were higher in samples collected in 2014 compared to those collected in 2015 on SCL (Supplementary Figure 5G). Samples collected from captive foxes had lower richness and evenness than those collected from wild foxes (Figure 4C).

Geographic source is a large contributor to bacterial community composition in U. littoralis scat (Figure 3). This result is expected, given that many biotic and abiotic factors change with biogeography which have been previously shown to affect the bacterial microbiome (Linnenbrink et al., 2013). Overall, the scat bacterial communities of U. littoralis were made up of the usual mammalian gut taxa; Bacteroidetes, Firmicutes, Fusobacteria, and Proteobacteria are consistent, abundant players in the canid scat bacterial communities (Ley et al., 2008; Zhang and Chen, 2010). In other studies, bacterial community differences among populations have been shown to be driven by differences in diet (Barelli et al., 2015), which could also contribute to the differences we see across islands in our dataset. For example, the lower relative abundance of Firmicutes classes Bacilli and Clostridia on SMI and SRI, and the higher abundance of Fusobacteriales, compared to the southern islands could reflect a more protein-rich diet (Deitch et al., 1987; Hang et al., 2012; Finlayson-Trick et al., 2017). This hypothesis is consistent with data showing that foxes on SRI consume the lowest amount of fruit among the islands and that on SMI the only fruit consumed by foxes is that of a non-native plant and instead they rely heavily on insects (Cypher et al., 2014). The month the sample was collected likely reflects differences in diet as well, as discussed later in this section. Interestingly, we found a unidirectional pattern in which a subset of samples from southern islands grouped with northern samples, but none of our predictor variables were able to explain this pattern (data not shown). More work is therefore needed to understand this finding and other deviations from a previous phylogenetic tree (Robinson et al., 2016). In future work we will therefore conduct a formal test of phylosymbiosis to better understand the impact of host genetics on the microbiome.

While geography had the greatest impact on the scat bacterial communities, weight, timing of the collection, and sex also contributed to the geographic patterns. Weight was a significant factor that influenced the bacterial community within and among islands, although as previously noted body size (including weight) is significantly different among islands (Wayne et al., 1991). Additionally, weight was a significant factor within an island, specifically on SCL, where the majority of the significantly differentially abundant ASVs in heavier samples were associated with Clostridia. This result is concordant with previous work that showed Clostridiales to be more abundant in obese pet dogs (Handl et al., 2013). Obesity has also been associated with increased levels of Firmicutes (Ley et al., 2005, 2006; Turnbaugh et al., 2009; Riva et al., 2017), and while foxes are rarely overweight especially in the wild, we found a similar pattern of increased amounts of Firmicutes in heavier foxes. Weight is a well-studied factor, especially in humans and mice, that affects the microbiome (Bäckhed et al., 2004; Ley et al., 2005, 2006; Riva et al., 2017), and thus also an expected result of our study.

The timing of sample collection also significantly influenced the scat bacterial communities in U. littoralis both among and within islands, specifically on SCL and SRI. The top NCBI Blast hit for the ASV with the largest effect (largest log2 fold change) in foxes in 2014 on SCL was Parasutterella secunda, which is proposed to be a core taxon that has a role in bile acid maintenance and cholesterol metabolism (Ju et al., 2019). The top NCBI Blast hit for the ASV with the largest effect in foxes in 2015 on SCL was Bacteroides coprocola, which is also associated with bile and the breakdown of plant polysaccharides (Kitahara et al., 2005). At the genus level Megamonas, shown to ferment glucose to short chain fatty acids (Chevrot et al., 2008; Sakon et al., 2008), consistently came up associated with collection timing. Temporal changes in the gut microbiome can be associated with seasonal changes in diet (Ren et al., 2017). The winter is typically the rainy season in southern California, which drastically changes the flora and fauna for the subsequent months. Exemplary of this process, the total rainfall on SRI and SCL in December of 2014 was 9–25 times greater than the surrounding months [Supplementary Figure 4; (Western Regional Climate Center, 2019)] and likely changed food source availability for the foxes. The top NCBI Blast hit for the ASV with the largest effect in foxes in November on SRI was Diplorickettsia massiliensis, which is an obligate bacterium of arthropods (Mediannikov et al., 2010). The top NCBI Blast hit for the ASV with the largest effect in foxes in December on SRI was Rosenbergiella australiborealis, which was isolated from plant nectar (Lenaerts et al., 2014). The taxa that correlate with collection timing, especially on SRI, appear to be more associated with environmental drivers of diet than with host metabolism.

Sex contributed to the differences in bacterial communities, which is consistent with several mammalian gut studies (McKenna et al., 2008; Phillips et al., 2012; Nelson et al., 2013; Menke et al., 2017). While this result was significant on two islands, the pattern was likely driven by one male sample on SNI so we focused on the differences found between sexes in CAT individuals. The top differentiated bacteria between sexes appear to be associated with diet. The top NCBI Blast hit for the ASV with the largest effect in females was Fusobacterium mortiferum, which metabolizes sucrose (Pikis et al., 2002). The top NCBI Blast hit for the ASV with the largest effect in males was Faecalibacterium prausnitzii, which is a well-studied gut commensal that produces butyrate and is indicated as a sentinel in gut health (Sokol et al., 2008; Khan et al., 2012). While sex-biased differences in diet can explain differences in gut microbiota (Bolnick et al., 2014), the previous study of the island fox’s diet did not look at differences between sexes (Cypher et al., 2014) so more work is needed here to tease apart that relationship. Alternatively, sex hormones have been proposed as a mechanism for sex-specific microbiome differences (Org et al., 2016).

Captivity significantly altered the scat microbiome of U. littoralis as evidenced by shifts in specific taxa (i.e., Firmicutes) and lower overall bacterial diversity. Despite the small sample sizes including biological replicates, this strong pattern emerged and is consistent with previous studies (Amato et al., 2013; Kohl et al., 2014; McKenzie et al., 2017; Guo et al., 2019; Schmidt et al., 2019). The top NCBI Blast hit for the ASV with the largest effect in wild foxes was Ruficoccus amylovorans, which is a diazotroph that reduces nitrate and degrades starch (Lin et al., 2017). In contrast, the top NCBI Blast hit for the ASV with the largest effect in captive foxes was Streptococcus alactolyticus, a potentially pathogenic species associated with diseases (ruminal acidosis, infective endocarditis, and colorectal cancer) in humans and animals (Jans et al., 2015; Jans and Boleij, 2018). Furthermore, several subspecies in the Streptococcus bovis/Streptococcus equinus Complex, to which S. alactolyticus belongs, have antimicrobial resistance (Jans and Boleij, 2018; Pompilio et al., 2019). While the mechanism of the bacterial shifts seen in this study remains unknown, previous studies have suggested changes in diet between captive individuals and their wild counterparts drive community differences in microbiota (Clayton et al., 2016; McKenzie et al., 2017; Metcalf et al., 2017; Greene et al., 2018; Oliveira et al., 2020). In the wild, island foxes are opportunistically omnivorous feeding on small mammals, insects, and native fruits (Cypher et al., 2014), while individuals in captivity are fed food items not found on the Channel Islands [e.g., dog food, hard-boiled eggs, fruits, and vegetables; (Coonan et al., 2010)]. Our results point to the influence that captivity can have on the fox immune system and subsequently the scat bacterial community. Reduced gut microbiome diversity can affect nutrient absorption and immune system function, ultimately affecting the overall health of the organism (Wienemann et al., 2011; Amato, 2013; Bragg et al., 2020). While there is some evidence that captive animals can gain back lost microbes upon release into the wild (Chong et al., 2019), carnivores have been shown to be more susceptible to losing microbial diversity in captivity (Metcalf et al., 2017) and there is evidence that it is more difficult to gain back wild microbes with an increasing number of generations in captivity (Webster et al., 2011; Becker et al., 2014). Going forward it will be imperative to manage host-associated microbiomes for captive foxes, especially if they will be released into the wild. To encourage a more diverse gut microbiome in captivity, managers could increase microbial reservoirs in enclosures (i.e., Loudon et al., 2014), shift diets to be more diverse and mimic natural food items, and decrease the use of antibiotics, as presented in Chong et al. (2019). More intensive approaches, including using prebiotics and probiotics that reflect the wild microbiome and/or fecal transplants could also be included to maintain diverse gut microbes.

For the captive samples, the PCoA (Figure 4A) and multivariate dispersion (Figure 4B) both show that SBZ2 and OCZ1 (both born in captivity) are similar to each other, while SCZ1 (born on SCL) is more similar to SCL. This pattern, while not confirmed statistically, suggests that inoculation from the natal environment may be more important than current environmental effects, which warrants further study. This legacy effect has been well-characterized in humans (Bokulich et al., 2016) and identified in other taxa (Jacob et al., 2015; Martínez-García et al., 2016; Gillingham et al., 2019). Some organisms have been able to maintain 60–70% of their wild microbes in captivity (Kohl and Dearing, 2014; Kohl et al., 2017), which could suggest that animals born in the wild but housed in captivity may have better health outcomes than those born in captivity.

In this first assessment of the scat microbiome across all islands of the Channel Island fox, we found patterns largely driven by geography, with food availability emerging as a strong contributing factor. These results add to the mounting evidence of distinctions among the six subspecies. Discussions of a genetic rescue to augment island populations of low genetic diversity by introducing animals from islands of higher genetic diversity will have to consider the microbiome differences as they represent possible divergences in natural history (diet) and immune systems that could be incompatible with other island-specific pathogens (Kohl et al., 2018). Secondly, even though there is extremely low genome-wide diversity in U. littoralis, a recent study found no signs of inbreeding depression in skeletal morphological traits (Robinson et al., 2018). The prevailing thought is that, on some of the islands, foxes have been able to persist despite the negative consequences of small population sizes and low genetic diversity. Perhaps other factors, such as the microbiome, contribute to their longevity and health, and the plasticity of the microbiome contributes to acclimatizing and eventually adapting to demographic and environmental changes (Alberdi et al., 2016). The insights provided in this study attest to the potential that microbiome research has for informing conservation science.