Microbiomes in Canidae

Abstract Because of their range expansion across North America, coyotes (Canis latrans) now occur sympatrically with numerous other predator species, including red foxes (Vulpes vulpes). This raises several interesting ecological questions, including if and how sympatry affects the diet and gut microbiomes of coyotes and red foxes. We examined the gut microbiomes of sympatric populations of coyotes and red foxes within two different National Parks in Virginia, USA, that differ in land use, vegetation, and anthropogenic disturbance: Prince William Forest Park (PRWI) and Manassas National Battlefield Park (MANA). From 2012 to 2017, scat samples from PRWI and MANA were collected and analyzed. Polymerase chain reaction (PCR) amplification of a region of the mitochondrial cytochrome‐b gene followed by restriction enzyme digestion of the PCR product was used to determine the origin of each scat sample. Next‐Generation DNA sequencing of a hypervariable 16S rRNA gene region was used to determine gut microbiome information about the scat samples. There was no evidence for a difference between the gut microbiomes of red foxes in either location, or for a difference between the gut microbiomes of red foxes at either location and coyotes at the location with lower human disturbance, PRWI. However, the gut microbiomes of coyotes at the location with higher anthropogenic disturbances, MANA, revealed a marked change from those found in red foxes at either location and from those in coyotes at the location with lower disturbances. The gut microbiomes of coyotes subjected to greater human impact may provide evidence of dysbiosis, indicative of increased physiological stress and reduced health. We discuss our observations in the context of understanding anthropogenic impacts on coyote and red fox interactions. Our results suggest that physiological stress in the form of human disturbance may play an important role in the composition of the gut microbiome of coyotes, which can affect their overall health.


| INTRODUC TI ON
There has been considerable interest in the ecological impact of the geographic expansion of coyotes (Canis latrans) across North America in both natural settings and anthropogenically altered systems. Because coyotes readily prey on domestic cats (Felis catus) (Crooks & Soule, 1999;Grubbs & Krausman, 2009), they can reduce the negative impact of feline predation pressure on vertebrate populations (Crooks & Soule, 1999). Considering that cats may kill up to 4 billion birds and 22.3 billion mammals in the United States annually (Loss et al., 2013), coyotes clearly could contribute to the maintenance of vertebrate diversity.
The presence of coyotes also has a predictable impact on red foxes (Vulpes vulpes). At the scale of the entire North American continent, interactions involving wolves, coyotes, and foxes result in coyotes being more abundant than foxes in areas without wolves, and foxes being more abundant than coyotes in areas where wolves are present, which are two examples of mesopredator release (Newsome & Ripple, 2015). In Minnesota, the reintroduction of wolves has had the effect of suppressing coyote populations to the benefit of fox populations (Levi & Wilmers, 2012). In areas where coyotes and foxes are sympatric in the absence of wolves, humans seem to be playing a mediating role, much as has been seen in coyote and cat interactions. Multiple studies have found that coyotes avoid human development, while foxes avoid coyotes and tend to concentrate in areas that are highly developed (Cove et al., 2012;Moll et al., 2018;Wait et al., 2018). Observations by Mueller et al. (2018) indicate that coyote/fox interactions are not necessarily antagonistic and harmful to foxes, but their data also show foxes and coyotes partitioning use of urban landscape, with foxes occupying more developed areas, and coyotes avoiding areas heavily used by humans.
Coyotes in poor health (where state of health was assessed by the visible presence of sarcoptic mange) were more likely to use human-developed areas (Murray et al., 2015), and the authors speculated that diseased coyotes might enter urban habitats in search of more easily obtained anthropogenic food sources.
Developments in DNA sequencing technology allow a unique and non-invasive approach to examining the health of animals in the wild, and the effects of environment and ecology through the analysis of their gut microbiomes (Bahrndorff et al., 2016;Barko et al., 2018;Ingala et al., 2018). Recent research has shown the microbiome to be critical to health in humans (The Integrative Human Microbiome Project Research Network Consortium, 2014), and this understanding has been extended to other species (Swanson et al., 2011). Increasingly, the gut microbiome is being studied to examine its role in health and disease in a range of non-human animals (Hanning & Diaz-Sanchez, 2015), providing new opportunities to explore the impact of a wide range of ecological interactions on wildlife. Here we describe our analysis of microbiomes of coyotes and foxes living in sympatry in two different National Parks representing very different levels of human habitation, development, and land use.

| Study sites
The study sites were Prince William Forest Park (PRWI) and Manassas National Battlefield Park (MANA), both located in Prince William County, Virginia, USA and part of the United States' National Park Service (NPS) system (Lookingbill et al., 2014). PRWI, an area of relatively low human disturbance, contains about 5075 hectares of undeveloped deciduous forest, with hiking trails and campgrounds present (Bozarth et al., 2011). MANA is about 30 km north of PRWI and comprises around 2100 hectares of open fields and forested areas (Rossell et al., 2005). A site of high human activity and disturbance, MANA contains many walking trails, open-air historic monuments, and picnic areas for visitors, has a major highway and other roads passing through it, and is surrounded by commercial and residential areas (Rossell et al., 2005).

| Scat collection
From 2012 to 2017, scat samples were collected across different seasons from both PRWI and MANA (permit PRWI-2013-SCI-0002).
Scat samples were collected either through the use of a grid system, which involved randomly choosing a plot area and then walking trails within that plot area to collect samples, or by walking various trails within the parks. Because aging of samples and potential contamination by microbes is difficult to assess in the field, we disregarded samples with visible microbial or fungi growth, and did not collect immediately after rain. A total of 264 scat samples were gathered.
Each scat sample was divided into subsamples for analyses in different labs at Towson University and stored at −20°C.

| DNA isolation
We used QIAamp DNA Stool Mini Kits (QIAGEN, Inc., Valencia, California, USA) to isolate total DNA from each scat sample. The Stool Pathogen Detection protocol, along with the recommended steps, was followed. Isolated DNA was stored at −20°C and used in downstream applications.

| Identifying canid species
To determine which canid species produced the scat samples, we followed a protocol described by Paxinos et al. (1997) to PCR amplify a ~400 bp region of the mitochondrial cytochrome-b gene.
Subsequent restriction enzyme digestion of the PCR product yielded an agarose gel banding pattern that was species-specific. To confirm the results of the restriction enzyme digestions, we used a protocol developed by Bozarth et al. (2010). We PCR amplified a portion of the mitochondrial d-loop region and the resulting PCR products were sized on agarose gels. In this case, PCR product size allowed species identification.

| Determining the gut microbiome
To determine the bacteria community present in the gut of the canid species, we used primers to PCR amplify a hypervariable V3 and V4 region of the 16S rRNA gene (Klindworth et al., 2013) in DNA isolated from the scat samples. We followed the protocol described in Illumina's 16S Metagenomic Sequencing Library Preparation guide (https://suppo rt.illum ina.com/downl oads/16s_metag enomic_seque ncing_libra ry_prepa ration.html).

| DNA sequencing
Following amplification of a region of the 16S rRNA gene for all samples, we used a unique pair of index primers for each sample in a second PCR step, following the protocol described in Illumina's 16S Metagenomic Sequencing Library Preparation guide. We then pooled samples and sequenced them on an Illumina MiSeq Next-Generation DNA sequencer using a 2x250 cycle Nextera XT sequencing kit (both from Illumina, Inc., San Diego, California, USA) and following the protocol described in Nextera XT Preparing DNA Libraries for Sequencing on the MiSeq (https://suppo rt.illum ina.com/downl oads/ prepa re_libra ries_for_seque ncing_miseq_15039 740.html).
Of note, DNA isolation and the first PCR set-up were both done in a different laboratory from all subsequent post-PCR activity to avoid potential contamination.

| Analyses
We used Illumina's 16S Metagenomics application found on BaseSpace to analyze our sequencing data. This application compares the sample sequences to reference sequences in the Greengenes database to classify the bacteria at different taxonomic levels (https:// suppo rt.illum ina.com/downl oads/16s_metag enomic_seque ncing_ libra ry_prepa ration.html). For each sample, this application reports the percent of sequence reads corresponding to a specific bacterial taxonomic group. In this study, the identification of the bacteria is reported at the phylum level and family level.
The sequencing data, specifically the percent of sequence reads in a sample that correspond to a specific bacterial phylum, were compared among the same canid species between the two different parks (PRWI and MANA) and between the two different canid species (coyote and red fox) within the same park. Multivariate analyses involving distance measures were used to determine if there was an association between the gut microbiome and a particular canid species and park, with a significant p-value of .05 or less (McCune & Grace, 2002;Suchodolski, 2016). We used the software program PC-ORD to carry out Multi-response Permutation Procedures (MRPP), which is a nonparametric procedure that can be used to test whether a significant difference exists between two groups of sampling units (the scat samples in this case) (Shin et al., 2015), allowing for the comparison of gut microbiome between canid species and park. PC-ORD also calculates the species richness (S), species evenness (E), Shannon's diversity index (H), and Simpson's diversity index (D) for various group (canid and park) comparisons. In addition, we used Kruskal-Wallis rank sum tests to compare the Shannon diversity and microbial species richness between host species and park.
Illumina's 16S Metagenomics application provided the Shannon diversity index for each sample, as well as the number of microbial species present in each sample; these values were used in the Kruskal-Wallis rank sum tests.
For the gut microbiome analyses using MRPP, two matrices were created, with the main matrix containing rows for the scat samples and columns for the bacterial phyla, which are categorical variables.
Seven bacterial phyla were included, which were Actinobacteria,

| Samples
The DNA isolations were performed on a total of 198 scat samples collected from PRWI and MANA. Of note, 66 of the 264 collected samples did not undergo the DNA isolation procedure; these samples in many cases contained an insufficient amount of digested material for the procedure (e.g., consisted of lots of hair, bone, or other material hard to digest). Of the 198 scat samples subjected to the DNA isolation procedure, a total of 161 samples were successfully PCR-amplified, sequenced, and subsequently analyzed to determine the gut microbiome of coyotes and red foxes from the two parks. Of note, 37 of these 198 scat samples could not be PCRamplified despite repeated attempts, possibly because the DNA isolation procedure with these samples was unsuccessful due to sub-optimal amounts of scat being used during the procedure or possibly because these samples contained degraded bacterial DNA.
Of the 161 scat samples that were analyzed, 58 were produced by coyotes, and 67 were from red foxes. Three scat samples were from gray foxes (Urocyon cinereoargenteus), eight samples came from domestic dogs (Canis lupus familiaris), one sample was from a domestic cat (Felis catus), and one sample was from a bobcat (Lynx rufus). The source of 23 of the 161 samples could not be identified.

| Coyote gut microbiome results
Of the 58 coyote gut microbiome samples, 37 were from PRWI and 21 from MANA. Figure 1 shows the percent of reads corresponding to each bacterial phylum in the coyote scat samples from each park

| Red fox gut microbiome results
Of the 67 red fox gut microbiome samples, 27 were from PRWI and 40 were from MANA. Figure 1 shows the percent of reads corresponding to each bacterial phylum in the red fox scat samples from each

| Comparison of canid gut microbiomes
Overall, when considering the samples individually, the most abun- indicating that the differences observed in gut microbiomes between coyotes and red foxes in MANA were unlikely to be due to chance.
In addition, based on Kruskal-Wallis rank sum tests, the average Shannon diversity of the gut microbiome was not demonstrably different between hosts or sites, with p-values > .05. Also, based on these tests, the number of microbial species in the scat samples did not differ significantly between hosts (p-value > .05), but it did differ significantly between sites, with a p-value << .05 (Tables 1 and 2).  Note: The distance measures and groups used, as well as the values for the average distance within group, the test statistic (T), the observed delta (δ), the expected delta (δ), and the chance-corrected within-group agreement (A), and the p-value are shown for each gut microbiome comparison. The gut microbiome of coyotes between PRWI and MANA, red foxes between PRWI and MANA, coyotes and red foxes in PRWI, and coyotes and red foxes in MANA was analyzed. The distance measures used were Euclidean (Pythagorean), Jaccard, and Sorensen (Bray-Curtis).

| DISCUSS ION
Based on our analyses of microbiome composition, there was no support for a difference between the gut microbiomes of red foxes at MANA and PRWI, nor between the gut microbiomes of red foxes at either location and coyotes at PRWI. By contrast, the composition of the gut microbiomes of coyotes at MANA were substantially different from the gut microbiomes of coyotes at PRWI and the red foxes at MANA and PRWI. Specifically, the MANA coyote samples were marked by a substantial increase in Proteobacteria, and a substantial decrease in Firmicutes. Firmicutes are typically found in the gut microbiota of humans, dogs, and other mammals, and a decrease in abundance has been proposed as an indicator of dysbiosis (Barko et al., 2018;Litvak et al., 2017). By contrast, while Proteobacteria are also commonly found in the guts of humans, dogs, and other mammals, an increase in Proteobacteria has been proposed as an indicator of dysbiosis (Blanchet et al., 2009;Hollins & Hodgson, 2019 (Gao et al., 2018).
We have observed that coyotes (but not red foxes) in areas more heavily utilized by humans show gut microbiome shifts associated with poor health. Multiple studies indicate that coyotes (but not red foxes) prefer to avoid heavily human-impacted areas (Cove et al., 2012;Moll et al., 2018;Mueller et al., 2018;Wait et al., 2018). It seems likely that coyotes forced into such areas are physiologically stressed, which can result in the gut microbiome changes found in multiple species (Hollins & Hodgson, 2019), which in turn can induce disease. Indeed, experimental manipulation of laboratory mice indicates that many of the deleterious effects of stress are actually caused by changes in the gut microbiome (Gao et al., 2018). Further, genetically identical mice raised under identical conditions, but with differing gut microbiomes, show differing resistance to infectious disease (Velazquez et al., 2019). It is therefore possible that what we have observed is the result of a concatenation of circumstances: coyotes forced into areas they consider sub-optimal are stressed, which leads to changes in the gut microbiome, which in turn is causal of a diseased state, including a damaged immune response and poor absorption of nutrients. This model is further supported by our data on red foxes. There is no reason to assume red foxes are stressed by occupying areas heavily utilized by humans; multiple studies found that red foxes seek out highly developed areas when coyotes are present (Cove et al., 2012;Moll et al., 2018;Mueller et al., 2018). seen in coyotes is caused by the presence of anthropogenic toxins, such as anticoagulant rodenticide that has been reported to affect urban bobcat populations (Riley et al., 2007;Serieys et al., 2018), and that has been suggested as a possible cause of anthropogenic stress in coyotes (Murray et al., 2015).
Increasingly, the gut microbiome is being shown to play a highly significant role in health and disease. Here we compared two competing species of canids with marked differences in toleration of human activity. In so doing, we provided a greater insight into the possible significance of the gut microbiome in the health of wild populations (Kartzinel et al., 2019). Our data suggest a working model that makes only a few assumptions that are supported by the literature. Coyotes forced into areas heavily utilized by humans become stressed, altering their gut microbiomes, resulting in poor health, poor nutrition, and an increased susceptibility to disease. We are not suggesting that this model can explain all the results observed when wild animals forced into urban habitats are studied, but we think that it seems likely that the gut microbiome may play a central role in a positive feedback cycle of stress and decreased resistance to disease that should be considered when studying any perturbed population.

This research was supported by a grant from the United States
National Park Service (P13AC00443). We thank Allan O'Connell for logistical support, as well as Katherine Morris and many volunteers for many hours in the field. We thank Towson University for providing a research grant and the facilities to perform the research. We thank David Hearn, Susan Gresens, and Nadim Alkharouf for assisting with the statistical analyses.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
We provide summary data within the manuscript text. Selected course material can be requested directly from the authors.