\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{The gut microbiomes of Channel Island foxes and island spotted skunks
exhibit fine-scale differentiation across host species and island
populations}



\author[1]{Samantha Pasciullo Boychuck}%
\author[2]{Lara Brenner}%
\author[3]{Calypso Gagorik}%
\author[4]{Juliann Schamel}%
\author[4]{Stacy Baker}%
\author[5]{Elton Tran}%
\author[5]{Bridgett von Holdt}%
\author[6]{Klaus-Peter Koepfli}%
\author[7]{Jesús Maldonado}%
\author[1]{Alexandra DeCandia}%
\affil[1]{Georgetown University}%
\affil[2]{The Nature Conservancy}%
\affil[3]{Northern Arizona University}%
\affil[4]{National Park Service}%
\affil[5]{Princeton University}%
\affil[6]{Smithsonian-Mason School of Conservation}%
\affil[7]{Smithsonian's National Zoo and Conservation Biology Institute}%


\vspace{-1em}



  
  \date{May 25, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
California's Channel Islands are home to two endemic mammalian
carnivores: island foxes (Urocyon littoralis) and island spotted skunks
(Spilogale gracilis amphiala). Although it is rare for two insular
terrestrial carnivores to coexist, these known competitors persist on
both Santa Cruz Island and Santa Rosa Island. We hypothesized that
examination of their gut microbial communities would provide insight
into the factors that enable this coexistence, as microbial symbionts
often reflect host evolutionary history and contemporary ecology. Using
rectal swabs collected from island foxes and island spotted skunks
sampled across both islands, we generated 16S rRNA amplicon sequencing
data to characterize their gut microbiomes. While island foxes and
island spotted skunks both harbored the core mammalian microbiome, host
species explained the largest proportion of variation in the dataset. We
further identified intraspecific variation between island populations,
with greater differentiation observed between more specialist island
spotted skunk populations compared to more generalist island fox
populations. This pattern may reflect differences in resource
utilization following fine-scale niche differentiation. It may further
reflect evolutionary differences regarding the timing of intraspecific
separation. Considered together, this study contributes to the growing
catalog of wildlife microbiome studies, with important implications for
understanding how eco-evolutionary processes enable the coexistence of
terrestrial carnivores -- and their microbiomes -- in island
environments.%
\end{abstract}%



\sloppy


\textbf{Introduction}

In comparison to mainland habitats, islands often harbor significantly
less vertebrate diversity . This is likely due to the difficulty of
island colonization, which can necessitate crossing large expanses of
water or being introduced by humans . The prevalence of insular
mammalian carnivorans (Order Carnivora) is especially low due to factors
such as large body size, low carrying capacity, reduced vagility, and
high extinction rate, with reduced space and limited prey availability
rendering niche differentiation difficult . It is therefore rare for two
endemic mammalian carnivores to coexist on the same island.

Counter to these ecological patterns, Channel Island foxes
(\emph{Urocyon littoralis santacruzae and U. l. santarosae} ) and island
spotted skunks (\emph{Spilogale gracilis amphiala} ) coexist on two
islands off the coast of southern California: Santa Cruz Island and
Santa Rosa Island. Both mesocarnivore taxa are endemic to the Channel
Islands, where they have coexisted alongside one another for at least
7,000 years . While island foxes are found on six of eight islands in
the archipelago, island spotted skunks are only found on the two largest
islands. Throughout their range, both species have undergone significant
population fluctuations and currently exhibit low levels of genetic
diversity, rendering them of particular conservation concern . These
species are also of interest to population ecologists, as they serve as
a rare example of mammalian carnivore coexistence in island habitats.

This coexistence occurs amid known competition between island foxes and
island spotted skunks. For example, direct competition was observed in
the early 1990s following the colonization of two invasive species to
the Channel Islands . As populations of feral pigs increased, they
served as abundant sources of prey that allowed golden eagles, an apex
predator from the mainland, to colonize the archipelago. Golden eagles
altered the trophic levels of the islands, turning the resident
predators -- namely, island foxes -- into prey . This intense predation
on island foxes changed the competitive hierarchy between the island's
endemic mesocarnivores by almost entirely removing island foxes from the
ecosystem. In response, island spotted skunk numbers increased and they
became more generalist in their diet and habitat use . The disruption of
competitive dynamics between island foxes and island spotted skunks
highlights the fragility of trophic structures within the Channel
Islands. Additionally, the decrease in island fox populations and
subsequent niche broadening of island spotted skunks demonstrates that
foxes are likely the dominant competitor and that fine-scale niche
partitioning enables coexistence between these two terrestrial
carnivores .

Island foxes and island spotted skunks partition resources in numerous
ways. While their niches overlap in habitat, temporal use, and diet,
they maintain subtle, but important, distinctions. Regarding habitat,
island foxes and island spotted skunks cohabitate within the coastal
Mediterranean ecosystem on Santa Rosa and Santa Cruz Islands. Yet,
island spotted skunks tend to be more specialized in their habitat use,
while island foxes utilize all areas across the island . There may be
further differences in microhabitat use, with island spotted skunks
often found in areas with ground-level cover in the form of rugged
terrain or dense vegetation . Regarding temporal dynamics, island
spotted skunks are largely nocturnal, whereas island foxes are
cathemeral (active during both the day and night; , with differing
degrees of seasonality. Finally, island spotted skunks are strictly
carnivorous while island foxes are largely omnivorous . Thus, island
spotted skunks tend to be more specialist and island foxes more
generalist along numerous ecological axes.

Another way to explore niche differentiation between island foxes and
island spotted skunks is to consider the composition of their
host-associated microbiomes. Microbiomes consist of all the internal and
external bacterial taxa that have coevolved alongside their host . These
bacteria function in modulating metabolism , responding to disease ,
aiding in digestion , mediating stress , and determining behavior . As
such, characterizing the composition of host-associated microbiomes can
provide an array of information about the host species. It can reflect
host diet , health status , environmental context , behavior , and even
evolutionary history ; Figure 1). While specific bacterial taxa often
differ between host species, the functional role of bacteria is largely
conserved among closely related hosts . For example, carnivores require
different bacteria than herbivores to assist in digestion, but both
ultimately help in the breakdown of food. Thus, the interplay between an
organism's phylogenetic history and external environment manifest in the
bacterial composition of their microbiome . Through understanding these
host-associated microbiomes, we can gain greater insights into the
functioning of the whole organism, their evolutionary history, and how
they fit within their broader ecosystem.

In the present study, we explored the role of gut microbial communities
in niche differentiation between island foxes and island spotted skunks
inhabiting Santa Cruz and Santa Rosa Islands. More specifically, we
examined the relative contributions of evolutionary history (species)
and environmental context (island) by comparing microbial variation
across host species and island populations using 16S rRNA amplicon
sequencing data. We predicted that the microbiomes of these two species
would exhibit high-level similarity due to their shared phylogenetic
ancestry as members of the suborder Caniformia within the mammalian
order Carnivora , low levels of genetic diversity, and overlapping
habitats across two islands. However, we also anticipated finer-scale
specialization of microbes at lower taxonomic levels, following
divergence of their ancestral lineages roughly 49 million years ago and
more recent niche differentiation through subtle differences in diet,
space use, and temporal activity patterns. Finally, we hypothesized that
comparatively specialist island spotted skunks would exhibit lower
microbial variation than more generalist island foxes, particularly in
light of skunks' smaller population sizes. Considered together, these
results provide important insights into the unusual coexistence between
these two insular mammalian carnivores. This not only broadens our
understanding of the eco-evolutionary dynamics operating between hosts
and their microbes in wild populations, but may further inform
conservation monitoring and management of at-risk species and the
symbionts they harbor.

\textbf{Materials and Methods}

\emph{Study Area}

The Channel Islands have a coastal Mediterranean climate with rainy
winters and dry summers and are made up of rocky coasts, grasslands,
shrublands, and woody forests . The islands have historically been
subject to disruption, principally by humans, through the introduction
of livestock and other invasive animal and plant species. This in turn
has heavily impacted native flora and fauna communities on the islands.
Santa Cruz Island (250 km\textsuperscript{2}, largest of the Channel
Islands) and Santa Rosa Island (215 km\textsuperscript{2}, second
largest) are part of the northern Channel Island archipelago and
separated by a 9.6 km wide channel. The two islands were once connected
as a super island called Santarosae that was separated from the mainland
until sea level rise separated them from one another approximately
10,000 years ago after the end of the Pleistocene epoch . Santa Cruz
Island (SCZ) is characterized by two east-west mountain ranges that
bound a central valley, and hosts a variety of habitat types including
grassland, chaparral, coastal sage scrub, oak woodland, and pine forests
. Santa Rosa Island (SRI) is characterized by a single mountain chain
running east-west through the middle of the island, with relatively flat
marine terraces bisected by canyons and dominated by grassland on the
north side, and more rugged terrain dominated by sagebrush scrub on the
south side. Small pockets of native island oak and pine woodland are
present. Both islands support a number of similar habitat types;
however, Santa Rosa Island is less rugged and less topographically and
vegetatively diverse compared to Santa Cruz Island. Grazing livestock
were removed from most of Santa Cruz Island by the late 1980s, while
grazers persisted on Santa Rosa Island until 2011. Thus, the majority of
Santa Cruz Island is further advanced in vegetative recovery than Santa
Rosa Island.

\emph{Sample Collection}

We sampled island foxes and island spotted skunks following live capture
in Tomahawk box traps during annual island-wide carnivore monitoring
surveys conducted between June 2021 and February 2022. We inserted a
Passive Integrated Transponder (PIT) tag in each individual animal for
identification purposes and recorded basic demographic information
including sex, weight, body condition, reproductive status, and
approximate age for each animal. We used sterile BD
BBL\textsuperscript{} CultureSwab\textsuperscript{}~swabs (BD
Biosciences, New Jersey, USA) to collect gut microbiome samples from the
perianal area by gentling swabbing a total of 40 times per animal
(\emph{i.e.} , 10 times before rotating the swab 90 degrees). All field
capture and handling efforts were conducted in accordance with
island-specific permissions and permits and following the guidelines
established by the American Society of Mammalogists for the ethical use
of wild mammals in research . Samples on Santa Cruz Island were
collected under a MOU with the California Department of Fish and
Wildlife (CDFW) and Scientific Collecting Permit (No. 008981) issued by
CDFW. Samples on Santa Rosa Island were collected under Channel Islands
National Park Research Permit \# CHIS-2023-SCI-0017. All protocols were
further approved by the Princeton University Institutional Animal Care
and Use Committees (Princeton IACUC \#3073).

\emph{DNA Extraction, PCR, and Sequencing}

We extracted microbial DNA from perianal swabs using a modified Qiagen
DNeasy PowerSoil HTP 96 Kit protocol (Qiagen Inc., Maryland, USA). We
isolated swab tips using scissors sterilized by dunking blades in 80\%
ethanol solution before passing them through a flame. We randomized
samples from different islands and species across PowerBead plates to
minimize potential batch effects, and disrupted all samples and controls
using a MiniBeadBeater-96. We shook plates for 3min, allowed 30sec of
rest, rotated plates, and shook plates for an additional 3min. We then
followed the standard manufacturer protocol before eluting samples for
15min at room temperature in 60\selectlanguage{greek}μ\selectlanguage{english}L solution C6 pre-heated to 70\selectlanguage{ngerman}ºC.

We used polymerase chain reaction (PCR) to amplify and tag the 16S rRNA
hypervariable 4 (V4) region with the barcoded forward
(GTGCCAGCMGCCGCGGTAA) and reverse (TAATCTWTGGGVHCATCAGG) primers
described in . PCR reactions contained 5\selectlanguage{greek}μ\selectlanguage{english}L HiFi HotStart ReadyMix (KAPA
Biosystems, Massachusetts, USA), 3.2\selectlanguage{greek}μ\selectlanguage{english}L primer mix (1.25\selectlanguage{greek}μ\selectlanguage{english}M), and 2.0\selectlanguage{greek}μ\selectlanguage{english}L
template DNA. Extraction controls included sterile swab tips (negative
controls) and DNA extracted from ZymoBIOMICS Microbial Community
Standard (D6300; positive controls); and PCR controls included PCR
reagents with no template DNA added (negative controls) and 1\selectlanguage{greek}μ\selectlanguage{english}L of
ZymoBIOMICS Microbial Community DNA Standard (D6305; positive controls;
Zymo Research, California, USA). PCR cycling conditions included:
initial denaturation of 94\selectlanguage{ngerman}ºC/3min, touchdown cycling for 30 cycles of
{[}94ºC/45s, 80-50ºC/60s, 72ºC/90s{]} decreasing 1ºC each cycle, 12
cycles of {[}94ºC/45s, 50ºC/60s, 72ºC/90s{]}, and a final extension of
72ºC/10min. To maximize yield and minimize PCR artifacts, each sample
was run through PCR in duplicate and subsequently combined for
quantification using a Qubit Flex Fluorometer (Invitrogen,
Massachusetts, USA). We then pooled 40ng of each barcoded library, used
the QIAquick PCR Purification kit (Qiagen Inc., Maryland, USA) and E-Gel
2\% Size Selection Protocol (Invitrogen, Massachusetts, USA) to clean
and size select libraries, and subsequently used Agencourt AMPure XP
magnetic beads (Beckman Coulter, California, USA) to concentrate final
libraries using a 1.25X bead to library ratio. Libraries were paired-end
sequenced (2x150nt) across three sequencing plates on an Illumina MiSeq
in the Princeton University Lewis Sigler Genomics Core Facility.

\emph{Sequence Processing and Analysis}

We used a paired-end, dual-indexed barcode splitter implemented
in\emph{Galaxy} to demultiplex raw sequencing data, with allowance for
one nucleotide mismatch between expected and observed barcodes. We then
imported demultiplexed data into QIIME 2 v2021.11 ; https://qiime2.org)
and ran each sequencing plate through the \emph{dada2
denoise-paired}plugin to further process reads. This pipeline corrected
likely sequencing errors, trimmed low quality bases, merged paired-end
reads, and removed chimeric reads. We then combined data from all three
sequencing plates into one dataset for further filtering based on reads
that appeared in positive and negative control samples.

We included six positive controls to assess possible sources of bias
during DNA extraction (n=3, ZymoBIOMICS Catalogue No. D6300) and PCR
amplification (n=3, ZymoBIOMICS Catalogue No. D6305). All six samples
contained sequences from the three Gram-negative and five Gram-positive
bacterial species contained within the mock community standards. This
confirmed successful lysis, DNA extraction, amplification, and taxonomic
identification of bacteria contained within our dataset. To identify
potential contaminants, we implemented the R package
\emph{decontam}using default parameters for the combined method that
considers the frequency (distribution of the frequency of each sequence
feature) and prevalence (presence/absence across samples) of potential
contaminants, with sequencing plate used as a batch identifier. We
subsequently filtered our dataset to remove the 32 amplicon sequence
variants (ASVs) identified as possible contaminants by \emph{decontam} ;
130 ASVs with taxonomic labels of ''mitochondria'' or ''chloroplasts'';
any ASVs with a total frequency \textless{}10; and any ASV found in only
one sample. We then removed all positive controls, negative controls,
duplicate samples, and samples with fewer than 2500 reads. This produced
a final dataset of 3,669,688 sequences and 1,443 ASVs sequenced in 110
samples, including 21 foxes and 22 skunks sampled on Santa Cruz Island,
and 33 foxes and 34 skunks sampled on Santa Rosa Island (Table S1).

\emph{Statistical Analyses}

We used QIIME 2 version 2022.8 ; https://qiime2.org) to perform all
downstream statistical analyses. We first implemented
the\emph{core-metrics-phylogenetic} function with sampling depth set to
2500 and a midpoint rooted phylogeny created with the
\emph{alignment}and \emph{phylogeny} functions in order to calculate
multiple measures of diversity. Alpha diversity metrics included
observed features to measure richness , Pielou's evenness to measure
species equitability , and the Shannon Diversity Index as a combined
measure of richness and evenness . Beta diversity metrics included
Bray-Curtis Dissimilarity to measure differences in microbial abundance
and the Jaccard Index to measure differences in species presence . We
calculated these five diversity measures in order to tease apart which
aspects of diversity (\emph{e.g.} , bacterial presence vs. abundance)
differed between our groups of interest. We performed multivariate
significance testing using the \emph{longitudinal anova} function for
alpha diversity measures and the \emph{diversity adonis} function
(multivariate PERMANOVA; for beta diversity measures. These analyses
contained the following variables: island, species, their interaction,
and sequencing plate. We visualized beta diversity differences using
principal coordinate analysis implemented through the EMPeror plug-in .
All figures were created using the \emph{qiime2R} and \emph{ggplot2} R
packages.

We assigned taxonomy to representative ASVs by using a Naïve Bayes
classifier trained on Greengenes 13\_8 reference sequences , as
described previously . We then characterized the taxonomic composition
of island fox and island spotted skunk gut microbiomes using
the\emph{taxa barplot} function and identified differentially abundant
taxa within and between species using analysis of composition of
microbiomes (ANCOM; as implemented through the \emph{ancom} function.
ANCOM overcomes statistical hurdles associated with compositional
datasets by calculating pairwise log-ratios between all taxa present
within sampling groups. It identifies which taxa consistently differ in
abundance between groups by calculating the test statistic, \emph{W,} as
the number of times the null hypothesis of ``no difference'' is violated
for a particular taxon. As numerous differences are expected to arise
between host species, we performed this analysis at the taxonomic level
of class.

\textbf{Results}

\emph{Island foxes exhibited higher alpha diversity than island spotted
skunks across both Santa Cruz and Santa Rosa Islands}

We found that bacterial richness and Shannon diversity significantly
differed by species (ANOVA: richness, \emph{F=} 45.179, \emph{df}
=1,\emph{p} \textless{}0.001; Shannon diversity, \emph{F=}
5.777,\emph{df} =1, \emph{p=} 0.018) and sequencing plate (ANOVA:
richness,\emph{F=} 87.182, \emph{df} =2, \emph{p} \textless{}0.001;
Shannon diversity, \emph{F=} 7.559, \emph{df} =2, \emph{p}
\textless{}0.001), with evenness only differing by sequencing plate
(ANOVA: evenness,\emph{F=} 16.709, \emph{df} =2, \emph{p}
\textless{}0.001; Figure 2; Figure S1). In contrast, island and its
interaction with species did not significantly predict alpha diversity
across all three measures (all\emph{p} \textgreater{}0.05; Table S2). We
further explored these patterns using pairwise \emph{t-} tests
implemented between species and sequencing plates (Figure S1; Table S3).
Between species, skunks exhibited significantly lower richness (pairwise
\emph{t-} test:\emph{t} =-3.872, \emph{adj.p} \textless{}0.001; Figure 2
A) and lower Shannon diversity approaching significance (pairwise
\emph{t-} test:\emph{t} = -1.696, \emph{adj.p} =0.093; Figure 2 B).
Between sequencing plates, plate CIF\_S08 exhibited lower richness and
Shannon diversity amid higher evenness compared to plates CIF\_S09 and
CIF\_S10 (all\emph{adj.p} \textless{} 0.05), which did not significantly
differ from each other (all \emph{adj.p} \textgreater{} 0.05; Figure S1;
Table S3).

\emph{Bacterial composition was significantly influenced by host
species, with subtler inter-island differences more pronounced in island
spotted skunks compared to island foxes}

Both Bray-Curtis Dissimilarity and the Jaccard Index were significantly
predicted by species, island, their interaction, and sequencing plate
(adonis: Table 1). Across both beta diversity measures, species
explained the largest proportion of the variance in the dataset, with
24.0\% for Bray-Curtis and 19.1\% for Jaccard. Principal coordinate
plots mirrored this result (Figure 3 A-B), with separation by species
clearly evident along PC1, which explained 24.70\% of the variation for
Bray-Curtis and 20.49\% of the variation for Jaccard. Further
examination Bray-Curtis Dissimilarity revealed separation between foxes
by an unknown variable along PC2 (which explained 6.388\% of the
variation) and between skunks by island along PC3 (which explained
6.188\% of the variation; Figure S2). Using Jaccard Distances, we
observed separation between skunks by island along PC2 (which explained
5.985\% of the variation) and between sequencing plates along PC3 (which
explained 4.176\% of the variation; Figure S2).

Given the strong influence of species on bacterial community structure,
we created subsets of the data by species to look for finer-scale
intraspecific differentiation between island populations. Within island
spotted skunks (Figure 3 C-D), bacterial abundance and presence were
significantly predicted by island (adonis: Bray-Curtis, \emph{F}
=10.936,\emph{R\textsuperscript{2}} =0.159, \emph{df} =1, \emph{p=}
0.001; Jaccard,\emph{F} =10.386, \emph{R\textsuperscript{2}} =0.146,
\emph{df} =1,\emph{p=} 0.001) and sequencing plate (adonis:
Bray-Curtis,\emph{F} =2.967, \emph{R\textsuperscript{2}} =0.086,
\emph{df} =2,\emph{p=} 0.001; Jaccard, \emph{F}
=4.460,\emph{R\textsuperscript{2}} =0.125, \emph{df} =2, \emph{p=}
0.001). We observed similar patterns in island foxes (Figure 3 E-F),
although island (adonis: Bray-Curtis, \emph{F}
=3.625,\emph{R\textsuperscript{2}} = 0.063, \emph{df} =1, \emph{p=}
0.001; Jaccard, \emph{F} =3.379, \emph{R\textsuperscript{2}} =0.056,
\emph{df} =1,\emph{p=} 0.001) and sequencing plate (adonis:
Bray-Curtis,\emph{F} =1.934, \emph{R\textsuperscript{2}} =0.067,
\emph{df} =2,\emph{p=} 0.002; Jaccard, \emph{F}
=3.254,\emph{R\textsuperscript{2}} =0.109, \emph{df} =2, \emph{p=}
0.001) explained lower proportions of variation.

\emph{Island foxes and island spotted skunks exhibited the core
mammalian microbiome, with the relative abundance of bacteria differing
between host species and (for skunks) island populations}

Gut microbial communities of island foxes and island spotted skunks
contained high proportions of Firmicutes (SCZ: Fox=33.9\%, Skunk=43.4\%;
SRI: Fox=20.8\%, Skunk=56.8\%) and Bacteroidetes (SCZ: Fox=26\%,
Skunk=24.5\%; SRI: Fox=46.7\%, Skunk=25.6\%). Dominant classes within
these phyla included Bacilli, Clostridia, and (for foxes)
Erysipelotrichi (Phylum: Firmicutes); and Bacteroidia (Phylum:
Bacteroidetes; Figure 4). Additional phyla with relative abundances
above 1\% included Actinobacteria (SCZ: Fox=3.6\%, Skunk=21\%; SRI:
Fox=7.7\%, Skunk=5.8\%), Proteobacteria (SCZ: Fox=20.6\%, Skunk=4.1\%;
SRI: Fox=13.3\%, Skunk=7.1\%), and Fusobacteria (SCZ: Fox=15.1\%,
Skunk=5.6\%; SRI: Fox=9.7\%, Skunk=4.3\%). Dominant classes within these
phyla included Actinobacteria and Coriobacteriia (Phylum:
Actinobacteria); Gammaproteobacteria, Epsilonproteobacteria, and
Deltaproteobacteria (Phylum: Proteobacteria); and Fusobacteriia (Phylum:
Fusobacteria; Figure 4). All additional phyla summed to less than 2\%
within each island population (SCZ: Fox=0.8\%, Skunk=1.4\%; SRI:
Fox=1.8\%, Skunk=0.4\%).

While the presence of taxa was largely consistent across species and
islands, their relative abundances often differed. We therefore used
ANCOM to identify significantly differentially abundant taxa within each
island and species at the class level (Table 2). Within-island
comparisons between island foxes and island spotted skunks returned four
unique classes. On Santa Cruz Island, island foxes harbored more
Gammaproteobacteria (W=54; 13.8\% vs. 0.8\%) and Deltaproteobacteria
(W=50; 0.4\% v. 0.009\%) than island spotted skunks, while island
spotted skunks harbored more Actinobacteria (W=52; 9.2\% vs. 1.9\%) than
island foxes. On Santa Rosa Island, island foxes again harbored more
Gammaproteobacteria (W=55; 9.9\% vs. 0.4\%) and Deltaproteobacteria
(W=57; 1\% v. 0.01\%) than island spotted skunks, whereas island spotted
skunks exhibited more Epsilonproteobacteria (W=54; 6.5\% vs. 2\%).
Similar analyses conducted within each species revealed island-specific
differences within island spotted skunks but not within island foxes
(Table 2). Island spotted skunks exhibited a higher percentage of
Coriobacteriia (W=54; 11.8\% vs. 0.08\%) on Santa Cruz Island, with no
classes significantly more abundant in Santa Rosa Island skunks.

\textbf{Discussion}

The present study is among the first to compare the gut microbiome of
coexisting insular mesocarnivores within the context of fine-scale niche
differentiation. More specifically, we examined gut microbial diversity
and composition between island foxes and island spotted skunks
inhabiting Santa Cruz Island and Santa Rosa Island off the coast of
southern California. We found that both island foxes and island spotted
skunks harbor the core mammalian microbiome, while also exhibiting
species-specific microbial communities. Within species, island of origin
significantly influenced gut microbial communities, with more pronounced
differences observed in island spotted skunks compared to island foxes.
This distinction may reflect the more specialized nature of island
spotted skunk ecology, in contrast to the more generalist nature of
island fox ecology. It may further reflect evolutionary processes
operating during a longer period of separation between island spotted
skunk populations, which diverged roughly 10,000 years ago, compared to
island foxes, which colonized the islands roughly 7,000 years ago .

Regarding their similarities, the shared phylogenetic history of island
foxes and island spotted skunks within the class Mammalia, order
Carnivora, and suborder Caniformia contributes to both host species
harboring phyla that comprise the core mammalian microbiome. These phyla
include Firmicutes, Actinobacteria, Bacteroidetes, and Proteobacteria,
with the additional phylum Fusobacteria often reported in mammalian
carnivores . Interestingly, we observed higher levels of Fusobacteria in
island foxes, which is consistent with other studies characterizing the
microbiome of species within Canidae . These results reflect the
importance of evolutionary history in shaping the core features of
mammalian microbiomes at higher taxonomic levels .

Yet amid these high-level similarities, we also observed significant
differentiation between island foxes and island spotted skunks. In fact,
species identity explained the largest percentage of variation in our
dataset by a significant margin, with 24\% of microbial abundance and
19.1\% of microbial presence attributed to host species. This result may
reflect the evolutionary divergence between ancestral lineages of the
Canidae and Mephitidae families roughly 49 million years ago, as well as
ecological differences between these species that likely recur
throughout their range. We posit that the reported patterns of diversity
likely reflect both of these eco-evolutionary processes, and in
particular, the breadth of each species' ecological niche.

For example, examination of alpha diversity revealed a wider range of
bacterial ASVs detected in island foxes compared to island spotted
skunks, likely resulting from their more generalist resource utilization
. Island foxes encounter a wider variety of food sources and
environments that may contribute to their higher bacterial richness, as
microbiomes are known to reflect an organism's diet and environmental
setting . Conversely, as island spotted skunks are more specialized in
resource use, they likely experience reduced exposure to these external
sources of microbial variation, thus contributing to their lower
bacterial richness. Their diet provides a likely candidate underlying
this result. Compared to more omnivorous island foxes, island spotted
skunk diets are mostly carnivorous and may therefore require fewer
microbes to digest their smaller range of food sources . Previous
research further showed that carnivore microbiomes tend to exhibit lower
richness than omnivore microbiomes, as meat proteins are easier to
digest than plant compounds . We therefore concluded that species
identity exerted a strong influence on gut microbial communities, likely
due to their dietary differences.

These interspecific differences were further reflected in the relative
abundance of bacterial taxa, as numerous classes significantly differed
between island foxes and island spotted skunks. On both Santa Cruz and
Santa Rosa Islands, island foxes harbored higher relative abundances of
Gammaproteobacteria and Deltaproteobacteria than their island spotted
skunks counterparts. Gammaproteobacteria have been linked to coastal
amphipod species similar to those preyed upon by island foxes . They
have also been associated with soil invertebrates such as earthworms ,
which were recently introduced to the Channel Islands . These results
potentially reflect a dietary preference towards these food sources
among island foxes. It may also reflect evolutionary history, as
Gammaproteobacteria has been isolated from other species within Canidae
. The other class consistently more abundant in island foxes,
Deltaproteobacteria, has been linked to marine environments , and may
further reflect island fox foraging behavior on beaches .

We also noted interspecific differences unique to each island. On Santa
Rosa Island, island spotted skunks exhibit higher abundances of
Epsilonproteobacteria. The genus largely driving this result
was\emph{Campylobacter} , which is often associated with enteritis in
humans that come into contact with soil, untreated water sources, or
animals . In island spotted skunks, this may derive from their local
environment or prey items, as \emph{Campylobacter} has previously been
isolated from wild mice . On Santa Cruz Island, Actinobacteria proved to
be more abundant in island spotted skunks. As Actinobacteria is one of
the most abundant bacterial taxa in soils , this may reflect habitat use
by island spotted skunks, which spend more time in ravines or covered
places while foraging for soil-dwelling prey . We therefore posit that
the increased presence of certain bacterial classes in island foxes and
island spotted skunks may derive from eco-evolutionary differences in
their population histories, diets, and preferred habitat types.

While species identity was the strongest predictor of bacterial
composition, we also observed intraspecific differences between island
populations. This pattern was more pronounced in island spotted skunks,
where 15.9\% of bacterial abundance and 14.6\% of species presence were
attributed to island of origin, compared to 6.3\% of bacterial abundance
and 5.6\% of species presence in island foxes. We again attribute this
near three-fold difference to two likely factors: the longer separation
time and more specialist ecology of island spotted skunks, as subtle
differences between island habitats may have a stronger effect on their
gut microbial communities. This contrasts with the more recent
colonization events and more generalist ecology of island foxes, which
harbor a greater diversity of microbes that may reflect their
interactions with numerous dietary items and microhabitats across both
islands. This likely has a homogenizing effect on their gut microbial
communities, although we do still find evidence of intraspecific
variation, as previously reported in island fox gut microbial
communities .

Differential abundance testing supported these overarching patterns, as
intraspecific comparisons only revealed one bacterial class
significantly differing between island spotted skunk populations amid no
significant differences between island fox populations. Within skunks,
the Santa Cruz Island population harbored significantly more
Coriobacteriia, which has been linked with the addition of crickets to
the diet of domestic dogs as it helps break down chitin in insect
exoskeletons . These results suggest that Santa Cruz Island spotted
skunks may be eating larger proportions of native arthropod species than
their Santa Rosa counterparts. This link between Coriobacteriia and
crickets is particularly compelling, as previous research reported
Jerusalem crickets (\emph{Stenopelmatus fuscus} ) as the most frequently
consumed arthropod by Santa Cruz Island spotted skunks . Considered
alongside interspecific results, this finding further underscores the
importance of dietary preference as a driver of gut microbial
composition within and between host populations and species.

In summary, the present study characterized the gut microbiomes of
island foxes and island spotted skunks across two of California's
Channel Islands: Santa Cruz Island and Santa Rosa Island. We identified
the core mammalian microbiome present in both species likely due to
their shared evolutionary history, and reported evidence of niche
differentiation likely due to differences in diet and habitat use
between these coexisting mesocarnivores. We found that host species was
the primary driver of microbial differentiation within our dataset, and
additionally reported intraspecific variation between island
populations. These intraspecific differences were more pronounced in the
longer separated and more specialist island spotted skunks compared to
the more recently colonized and more generalist island foxes. These
results therefore support the claim that host evolutionary history and
contemporary ecology (\emph{i.e.} , diet, habitat use, and local
environmental context) all play a role in shaping the gut microbiome.
Commensal bacteria likely reflect -- and perhaps even contribute to --
the fine-scale niche differentiation that allows these insular
terrestrial carnivores to coexist. To better understand this ecological
rarity, we recommend further study of the diet, behavior, and bacteria
colonizing these known competitors. This will shed additional light on
the eco-evolutionary processes that have enabled their coexistence since
initial colonization millennia ago, and will hopefully facilitate future
coexistence for many millennia to come.

\textbf{Acknowledgments:} We would like to thank Bridget Parrino,
Heather Alford, Victor Zhang, and Mason London for collecting rectal
swabs from island foxes and island spotted skunks. We thank Robert
Fleischer and Nancy McInerney from the Center for Conservation Genomics
at Smithsonian's National Zoo and Conservation Biology Institute (NZCBI)
for their support of this project. This work was funded by the
Smithsonian Institution Fellowship Program, Friends of the Island Fox
Research Grant, Georgetown University Annual Research Grant, and
Georgetown University Summer Salary Supplement. Further support was
provided by the NZCBI Center for Conservation Genomics and Center for
Species Survival.

\textbf{References}

\textbf{Data Accessibility}

Sequencing data analyzed in this study is publicly available through the
NCBI Sequence Read Archive (SRA) under BioProject PRJNA1012705.
Demultiplexed forward and reverse reads are available for each sample
(BioSamples SAMN37271150 through SAMN37271259), with SRA accession
numbers SRR25897939 through SRR25898048. Metadata are also stored in the
SRA (BioProject PRJNA1012705) using the MIMARKS survey host-associated
6.0 package.

\textbf{Benefits Generated:} Benefits from this research accrue from the
sharing of our data and results on public databases as described above.

\textbf{Author Contributions}

ALD, LJB, CNG, JTS, SB, BMvH, KPK, and JEM designed the study. LJB, CNG,
JTS, and SB collected the samples. ALD conducted the laboratory work.
SPB, ET, and ALD processed and analyzed the data. SPB, LJB, CNG, JTS,
SB, ET, BMvH, KPK, JEM, and ALD prepared, contributed to, and approved
the final manuscript.

\textbf{Tables}

\textbf{Table 1.} Results from multivariate beta diversity analyses for
bacterial abundance (Bray-Curtis Dissimilarity) and presence (Jaccard
Index) as implemented through the \emph{diversity adonis} function
in\emph{QIIME2} . The \emph{F-} value (\emph{F} ), proportion of
variation explained (\emph{R\textsuperscript{2}} ), degrees of freedom
(\emph{df} ), and \emph{p} -value (\emph{p} ) are provided. Asterisks
indicate statistical significance.\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
& & F & R\textsuperscript{2} & df & p\tabularnewline
\midrule
\endhead
\textbf{Bray-Curtis} & \textbf{Island} & 6.420 & 0.040 & 1 &
0.001*\tabularnewline
& \textbf{Species} & 38.417 & 0.240 & 1 & 0.001*\tabularnewline
& \textbf{Plate} & 2.474 & 0.031 & 2 & 0.004*\tabularnewline
& \textbf{Island:Species} & 6.177 & 0.039 & 1 & 0.001*\tabularnewline
& \textbf{Residuals} & - & 0.650 & 104 & -\tabularnewline
\textbf{Jaccard} & \textbf{Island} & 6.030 & 0.040 & 1 &
0.001*\tabularnewline
& \textbf{Species} & 29.114 & 0.191 & 1 & 0.001*\tabularnewline
& \textbf{Plate} & 3.845 & 0.050 & 2 & 0.001*\tabularnewline
& \textbf{Island:Species} & 5.771 & 0.038 & 1 & 0.001*\tabularnewline
& \textbf{Residuals} & - & 0.682 & 104 & -\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 2. Five unique classes were significantly differentially
abundant between groups of interest.} ANCOM was conducted at the class
level within four subsets of data: Santa Cruz Island (SCZ) samples only,
Santa Rosa Island (SRI) samples only, island spotted skunks only, and
island foxes only. Across species comparisons returned four unique
classes, with Deltaproteobacteria and Gammaproteobacteria consistently
more abundant in foxes across both islands. Within species comparisons
returned one differentially abundant class between skunk populations,
with no significant results between island fox populations. The table
includes the Phylum and Class of each significant result alongside the
ANCOM test statistic (W, which indicates the number of times the null
hypothesis was rejected for a given class), the dataset used for the
analysis, and the direction of change (\emph{i.e.} , which group
exhibited higher abundance).\selectlanguage{english}
\begin{longtable}[]{@{}lllll@{}}
\toprule
\textbf{Phylum} & \textbf{Class} & \textbf{W} & \textbf{Dataset} &
\textbf{Direction}\tabularnewline
\midrule
\endhead
Actinobacteria & Actinobacteria & 52 & SCZ only & Foxes \textless{}
Skunks\tabularnewline
Proteobacteria & Deltaproteobacteria & 50 & SCZ only & Foxes
\textgreater{} Skunks\tabularnewline
Proteobacteria & Gammaproteobacteria & 54 & SCZ only & Foxes
\textgreater{} Skunks\tabularnewline
Proteobacteria & Epsilonproteobacteria & 54 & SRI only & Foxes
\textless{} Skunks\tabularnewline
Proteobacteria & Deltaproteobacteria & 57 & SRI only & Foxes
\textgreater{} Skunks\tabularnewline
Proteobacteria & Gammaproteobacteria & 55 & SRI only & Foxes
\textgreater{} Skunks\tabularnewline
Actinobacteria & Coriobacteriia & 54 & Skunks only & SCZ \textgreater{}
SRI\tabularnewline
- & - & - & Foxes only & SCZ = SRI\tabularnewline
\bottomrule
\end{longtable}

\textbf{Figure Captions}

\textbf{Figure 1. Fine-scale niche differentiation between island foxes
and island spotted skunks may influence their host-associated
microbiomes.} Although they exhibit shared evolutionary history (both
mesocarnivore species are within Order Carnivora) and inhabit the same
two Channel Islands (SCZ = Santa Cruz Island; SRI = Santa Rosa Island),
island foxes and island spotted skunks appear to coexist through niche
partitioning. Island foxes tend to be more generalist in their habitat
use, diet, and temporal activity, whereas island spotted skunks tend to
be more specialist in these areas. Commensal microbial communities may
reflect these differences, as depicted in this schematic created with
BioRender.com.

\textbf{Figure 2. Alpha diversity differs by species, but not
island.}(A) Bacterial species richness (as measured by observed
features) and (B) Shannon diversity significantly differed by species,
with foxes exhibiting higher diversity than skunks, but not by island.
(C) Bacterial evenness differed by neither species nor island.

\textbf{Figure 3. Beta diversity differs by species, island, and
sequencing plate.} We observed significant differences between species
when comparing (A) bacterial abundance (as measured using Bray-Curtis
Dissimilarity) and (B) bacterial presence (as measured using the Jaccard
Index). Differences between island populations (SCZ = Santa Cruz Island;
SRI = Santa Rosa Island) were greater within island spotted skunks (C,D)
compared to island foxes (E,F), with additional clustering evident by
sequencing plate. Fox and skunk silhouettes were created with
BioRender.com.

\textbf{Figure 4. Island foxes and island spotted skunks exhibited the
core mammalian microbiome, with differences in relative abundance
observed between species and island populations.} Each bar represents
the combined samples collected from each island population studied,
where SCZ indicates Santa Cruz Island and SRI indicates Santa Rosa
Island. Each color denotes a distinct bacterial class, with all
additional classes summed together and termed ``Other''. Taxonomic
groups annotated with an asterisk significantly differed between groups
of interest during differential abundance testing with ANCOM.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure1/Figure1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure2/Figure2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure3/Figure3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure4/Figure4}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Cover-Image-Bridget-Parrino/Cover-Image-Bridget-Parrino}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Cover-Image-John-Stuelpnagel/Cover-Image-John-Stuelpnagel}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

