Social environment and genetics underlie body site‐specific microbiomes of Yellowstone National Park gray wolves (Canis lupus)

Abstract The host‐associated microbiome is an important player in the ecology and evolution of species. Despite growing interest in the medical, veterinary, and conservation communities, there remain numerous questions about the primary factors underlying microbiota, particularly in wildlife. We bridged this knowledge gap by leveraging microbial, genetic, and observational data collected in a wild, pedigreed population of gray wolves (Canis lupus) inhabiting Yellowstone National Park. We characterized body site‐specific microbes across six haired and mucosal body sites (and two fecal samples) using 16S rRNA amplicon sequencing. At the phylum level, we found that the microbiome of gray wolves primarily consists of Actinobacteria, Bacteroidetes, Firmicutes, Fusobacteria, and Proteobacteria, consistent with previous studies within Mammalia and Canidae. At the genus level, we documented body site‐specific microbiota with functions relevant to microenvironment and local physiological processes. We additionally employed observational and RAD sequencing data to examine genetic, demographic, and environmental correlates of skin and gut microbiota. We surveyed individuals across several levels of pedigree relationships, generations, and social groups, and found that social environment (i.e., pack) and genetic relatedness were two primary factors associated with microbial community composition to differing degrees between body sites. We additionally reported body condition and coat color as secondary factors underlying gut and skin microbiomes, respectively. We concluded that gray wolf microbiota resemble similar host species, differ between body sites, and are shaped by numerous endogenous and exogenous factors. These results provide baseline information for this long‐term study population and yield important insights into the evolutionary history, ecology, and conservation of wild wolves and their associated microbes.

Despite increased study, there remain numerous questions about the primary factors underlying microbial species presence and abundance. Of particular interest is the contribution of host genetics, demography, and environment (Bonder et al., 2016;Ceja-Navarro et al., 2015;Goodrich et al., 2014;Kurilshikov et al., 2017;Rothschild et al., 2018;Spor et al., 2011). As microbiota function in diverse physiological processes, elucidating these factors can have important implications for the evolutionary history, ecology, and conservation of species (DeCandia et al., 2018;Hauffe & Barelli, 2019;Trevelline et al., 2019).
Within wildlife systems, studies often consider broad-scale patterns of phylosymbiosis or the eco-evolutionary scenario where host phylogenetic relationships are mirrored by dissimilarity between host-associated microbiomes (Brooks et al., 2016). This pattern has been observed within numerous host lineages, including invertebrates, rodents, and primates (Brooks et al., 2016). However, the taxonomic scale of comparison can strongly influence the degree of congruence between host phylogenies and microbial dendrograms, as these patterns break down within speciose host genera (Greene et al., 2019;Grond et al., 2020). Further, environmental and behavioral variables often explain significant portions of variance alongside phylogenetics. Metagenomic analyses of gut microbiomes across Mammalia, Aves, Reptilia, Osteichthyes, and others found evidence that host taxonomy, diet, lifespan, and behavior (i.e., activity and social structure) influenced microbial composition to varying degrees (Levin et al., 2021). Similar analyses of 16S rRNA across Mammalia, Aves, Reptilia, Amphibia, and Actinopterygii revealed that diet primarily predicted functional guilds and host phylogeny predicted the specific microbes present (Youngblut et al., 2019). The same pattern emerged within Mammalia, where diet was predictive of gut microbiome convergence at higher taxonomic rankings (such as microbial phylum) and host phylogeny was predictive of gut and skin microbial communities at lower taxonomic rankings (such as microbial family; Nishida & Ochman, 2018;Ross et al., 2018).
These broad-scale studies provide valuable insights into the evolutionary history of hosts and their associated microbes. However, they lack details that may be relevant to the ecology and conservation of lower host clades. In order to obtain finer-scale information, researchers have turned to species-specific studies. Within captive management settings, microbiome analyses have yielded important information about the reproductive (Southern white rhinoceros, Ceratotherium simum simum; Williams et al., 2019) and gastrointestinal (red wolf, Canis rufus; Bragg et al., 2020) health of captive-housed wildlife. Across host species, artificial diet and housing conditions significantly influence microbiota, as seen in primate microbes "humanized" by captivity (red-shanked douc, Pygathrix nemaeus, and mantled howler monkey, Alouatta palliata; Clayton et al., 2016 Gomez et al., 2015), with increasing interest in the contribution of within-species host genetics (house mouse, Mus musculus ;Suzuki et al., 2019).
The majority of in situ microbiome studies focus on gut microbiota (Gomez et al., 2015;Ingala et al., 2019;Sanders et al., 2015;Suzuki et al., 2019), with a few studies targeting skin sites (Avena et al., 2016;Hooper et al., 2019). While informative, narrow focus on one or two sample types provides an incomplete characterization of commensal microbes. Body sites exhibit distinct microbial communities, due to differences in oxygen exposure, nutrient availability, substrate, and environmental factors (The Human Microbiome Project Consortium, 2012). For example, human guts harbor far more anaerobic species than human skin, due to low oxygen availability in the gastrointestinal tract (Coates et al., 2019). Further differences exist between skin microenvironments, where dry sites (such as the forearm) exhibit higher species richness than sebaceous sites (such as the forehead; Byrd et al., 2018;Grice & Segre, 2011 Lavrinienko et al., 2018). This necessitates study of body site-specific microbiota in the wild to more holistically characterize host-associated microbiomes.
The present study bridges these knowledge gaps by characterizing in situ wildlife microbiomes across multiple body sites.
Specifically, we sequenced host-associated bacterial communities across six haired and mucosal body sites (and two fecal samples) in a wild pedigreed population of gray wolves (Canis lupus) inhabiting Yellowstone National Park (YNP) in Wyoming, USA (Figure 1). We hypothesized that body sites would harbor distinct microbial communities, with haired sites exhibiting higher alpha diversity than mucosal sites, and fecal sites hosting higher proportions of anaerobic bacteria than nonfecal sites. We further leveraged observational data, biobanked blood and tissue, and the highly resolved YNP wolf pedigree  to examine genetic, demographic, and environmental factors underlying body site-specific microbiota in this wild mammalian system. Given that wolves live in kin-structured family groups comprised of close relatives vonHoldt et al., 2008), we predicted that genetic relatedness would positively correlate with microbial similarity, as seen in studies of humans and model systems (Bonder et al., 2016;Spor et al., 2011). However, as many wolves disperse from their natal groups, we also predicted that unrelated wolves sharing the same social environment (here termed "sampling pack") would exhibit more similar microbiota than unrelated wolves sampled in different packs, as seen in cohabitating humans (Dill-McFarland et al., 2019;Song et al., 2013). Pack mates cooperatively share social contacts, resources, and prey items, including elk and (to a more variable degree) bison, deer, moose, or beaver -it therefore follows that they share microbes, as well. The social structure of wolves, coupled with the availability of host genetic, microbial, demographic, environmental, and pedigree information provided the unique opportunity to examine body site-specific microbiota and their underlying factors in a wild mammalian system. Results yielded system-specific insights, while also contributing to the larger-scale effort of characterizing wildlife microbiomes in situ.

| Sample and data collection
Gray wolves have been monitored annually by the National Park Service (NPS) since their reintroduction to YNP in 1995 and 1996.
Static (e.g., sex) and dynamic (e.g., pack membership) life history variables were collected from aerial and ground surveys conducted during the winter monitoring season. In addition, biological samples were collected between December and February during helicopter captures (whole blood sampled through venipuncture) and field necropsies (tissue sampled from deceased individuals). At the time of microbiome sample collection, body condition was qualitatively scored by two handlers based on sex-and age-specific patterns of weight, muscle/fat condition, coat condition, and injuries/illness detected. All capture and handling protocols were conducted in accordance with the NPS (IACUC permit IMR_YELL_Smith_wolves_2012) and Princeton University (Princeton IACUC #2009A-17) Institutional Animal Care and Use Committees.
We generated genomic data using whole blood and microbiome data using skin swabs and fresh scat. We used sterile BD BBL™ CultureSwab™ swabs to sample commensal bacteria at six body sites, including ear canal, nostril, lip commissure, axilla, dorsal flank, and perianal area ( Figure S1). At each body site, we rubbed the swab tip on the skin roughly 100 times, rotating by 90° every 25 times.
We also collected fresh scat if the individual released feces during or immediately prior to sampling. Across three field seasons, we col-

| Genomic DNA extraction, RAD sequencing, and data analysis
We extracted DNA from whole blood following the Qiagen DNeasy Blood and Tissue Kit manufacturer protocol (Qiagen Inc.), quantified DNA using Quant-iT™ PicoGreen™ dsDNA assays or high-sensitivity Qubit™ fluorometry (Thermo Fisher Scientific), visualized extracts on 1% agarose gels to assess molecular weight, and standardized concentrations to 5 ng/μl. We then generated genomic data using a modified restriction site-associated DNA sequencing protocol designed by (Ali et al., 2016) and described in (DeCandia et al., 2021;vonHoldt et al., 2020). Briefly, we digested DNA with the restriction enzyme sbfI and ligated uniquely barcoded, biotinylated adaptors. We pooled barcoded samples, sheared DNA to 400 bp using a Covaris LE220, and enriched for fragments containing the F I G U R E 1 Gray wolves were reintroduced to Yellowstone National Park in 1995-1996 and have been closely monitored ever since.

P H O T O
Credit: NPS/Jacob W. Frank ligated adaptor using a streptavidin bead-binding assay (Invitrogen Dynabeads M-280). We prepared libraries for sequencing using the NEBNext Ultra II DNA Library Preparation Kit (New England Biolabs) manufacturer protocol and performed size selection for fragments 300-400 bp in length using Agencourt AMPure XP magnetic beads (Beckman Coulter). We standardized final libraries to 10 nM and performed paired-end sequencing (2 × 150nt) on an Illumina HiSeq2500 or NovaSeq6000 at the Princeton University Lewis Sigler Genomics Core Facility.
After sequencing, we aligned all forward and reverse reads with the restriction enzyme cut site using a custom perl script (Data S1).
We then implemented the gstacks and populations modules in STACKS v2.2 (Rochette et al., 2019) to genotype and filter genomewide single nucleotide polymorphisms (SNPs) from paired-end data using the Marukilow model (Maruki & Lynch, 2017). We ran gstacks using a dataset of 32 samples representing 24 unique individuals (N.B., we excluded wolf 949M from these analyses due to microbial sampling after death). We manually determined which duplicate samples exhibited higher read counts and implemented populations using 24 unique samples and the filtering parameters --write_sin-gle_snp (which only retains one SNP per read) and -r 0.9 (which only retains loci genotyped in >90% of wolves). This yielded a dataset of 116,953 variant sites genotyped in 24 wolves. We subsequently removed singletons, doubletons, and X-chromosome sites (due to our mixed-sex sample set; Clayton, 2009) using VCFtools v0.1.12b (Danecek et al., 2011). This produced a dataset of 86,545 highconfidence autosomal SNPs found throughout the genome.
To consider familial relationships, we calculated pairwise relatedness coefficients in the R package related (Pew et al., 2015). We used the coancestry function and implemented the dyadic likelihood estimator (dyadml = 1; Milligan, 2003) with allowance for inbreeding (allow.inbreeding = TRUE). We additionally identified putative parents and grandparents of wolves with microbiome data by referencing the full YNP wolf pedigree, which included 871 parent-offspring pairs as of October 2020 (vonHoldt et al., 2020).

| Microbial DNA extraction and amplicon sequencing
We used a modified Qiagen DNeasy PowerSoil Kit protocol (Qiagen Inc.) to extract DNA from each sample, as described in (DeCandia et al., 2019. Briefly, we placed swab tips or fecal material into PowerBead tubes that were shaken for two cycles on a Qiagen TissueLyserII. Each cycle lasted for 12 min at 20 shakes/s, with the addition of 60 μl of C1 solution occurring between cycles. We then followed the standard manufacturer protocol until the final elution, when we incubated samples for 10-15 min at room temperature with 60 μl C6 buffer preheated to 70°C. We used sterile swab tips as negative controls during extraction and subsequent library preparation and sequencing. We concentrated extracts to 20 μl in a Vacufuge if needed, quantified DNA using a high-sensitivity Qubit™ fluorometer, standardized high-yield samples to 2.5 ng/μl, and included low-yield samples with concentrations as low as 0.062 ng/μl. We used barcoded forward (GTGCCAGCMGCCGCGGTAA) and reverse (TAATCTWTGGGVHCATCAGG) primers described in (Caporaso et al., 2011)  ear canal, n = 21; lip commissure, n = 24; nostril, n = 21; perianal area, n = 25; feces, n = 2; negative controls, n = 5) collected from 25 wolves (1-11 samples/wolf, with median = 5 and mode = 6). We used a barcode splitter for paired-end, dual-indexed data in the online platform Galaxy (Afgan et al., 2018) to demultiplex raw sequencing reads with allowance for a single nucleotide mismatch in sequence tags. We filtered 4,998,642 demultiplexed reads in QIIME 2 v2020.8 (Bolyen et al., 2019; https://qiime2.org) using the dada2 denoise-paired plugin (Callahan et al., 2016), which corrects probable sequencing errors, trims low-quality bases, merges paired-end reads, and removes chimeric sequences. We retained 4,272,465 sequences containing 8,944 amplicon sequence variants (ASVs ;   Table S1).

| Characterizing the host-associated microbiome across multiple body sites
We characterized microbial communities inhabiting each body site using QIIME 2 (Bolyen et al., 2019; https://qiime2.org). For these analyses, we employed four measures of bacterial diversity to consider different aspects of microbial community composition . For alpha (or within sample) diversity, we measured microbial species richness using observed amplicon sequence variants (ASVs; Hagerty et al., 2020) and species equitability using Pielou's evenness (Pielou, 1966). For beta (or between sample) diversity, we used qualitative (unweighted UniFrac) and quantitative (Bray-Curtis dissimilarity) measures to consider differences in species presence and abundance (Lozupone et al., 2007). Unweighted UniFrac is a qualitative measure that calculates the amount of branch length in a phylogenetic tree of ASVs leading to unique members of each microbial community (Lozupone & Knight, 2005). This measure primarily considers species presence and ignores the relative abundance of ASVs. In contrast, Bray-Curtis dissimilarity is a quantitative measure that directly incorporates abundance into its calculations (Bray & Curtis, 1957). By concurrently examining observed ASVs, Pielou's evenness, unweighted UniFrac distances, and Bray-Curtis dissimilarity, we were able to compare microbial species richness, evenness, presence, and abundance across body sites.
We calculated these measures using the core-metrics-phylogenetic and alpha-rarefaction functions in QIIME 2. We employed rarefaction to control for different sequencing depths and calculated unweighted UniFrac distances using a midpoint-rooted phylogeny constructed using the alignment and phylogeny functions. We used Kruskal-Wallis tests (Kruskal & Wallis, 1952) implemented through the alpha-group-significance function to compare observed ASVs and Pielou's evenness across body sites. To assess differences in beta diversity, we implemented PERMANOVA (Anderson, 2001) through the diversity adonis function (Oksanen et al., 2019) and performed principal coordinate analysis (PCoA) using the EMPeror plugin in QIIME 2 (Vázquez-Baeza et al., 2013).
In order to assign taxonomy to ASVs, we used reference sequences from the Greengenes 13_8 database to train a Naive Bayes classifier (Bokulich et al., 2018;DeSantis et al., 2006). Reference sequences were trimmed to match our target amplicon sequence, clustered at 99% similarity, and compared against representative sequences in our dataset using the classify-sklearn function in the feature-classifier plugin for QIIME 2. We subsequently clustered body sites hierarchically using Euclidean distances and the average linkage method. We reported the most abundant orders colonizing each body site to characterize high-level taxonomy and visualize broad-scale patterns of taxonomic diversity. We additionally considered lower-level taxonomy by implementing the linear discriminant analysis (LDA) effect size (LEfSe) method in Galaxy with sequences grouped at the genus level (Segata et al., 2011). This method identifies microbial clades (e.g., genera) that underlie differences between two or more biological groups (e.g., body sites). We used default parameter settings and the one-against-all strategy for multigroup analysis. Following LEfSe guidelines, we excluded the two fecal samples due to small within-group sample size (Segata et al., 2011).

| Assessing environmental, demographic, and genetic drivers of microbiome composition
We created two subsets of the microbiome data to perform body site-specific analyses. The first subset included perianal samples (n = 22) to represent the gut microbiome (Bassis et al., 2017) and the second subset included dorsal flank samples (n = 17) to represent the skin microbiome (Kong et al., 2017). We then considered potential environmental, demographic, and genetic drivers of microbial community composition. We tested for significant associations with age class, body condition, coat color, field season, sampling pack, and sex using single-and multifactor PERMANOVA implemented through the diversity adonis function in QIIME 2 (Anderson, 2001;Oksanen et al., 2019). We note that the variable field season was confounded with sequencing plate; thus, any observed differences may represent true signal of temporal change or artificial signal of technical batch effect.

TA B L E 1 Samples included in the microbiome dataset
We additionally tested for correlations between host genetic distance and microbial dissimilarity using full and partial Mantel tests implemented in the R package vegan with 999 permutations (Oksanen et al., 2019). Host genetic distances (Euclidean) were calculated using the dist function in the R package adegenet (Jombart, 2008;Jombart & Ahmed, 2011), and microbial distances were calculated in QIIME 2 using Bray-Curtis dissimilarity and unweighted UniFrac distances. We additionally constructed a dissimilarity matrix for pack membership (where 0 = same pack, 1 = different pack) to implement partial Mantel tests that controlled for shared sampling pack. The test statistic r was based on Spearman's rank correlation rho, and statistical significance was assessed using a threshold of p < .05.

| Genomic analyses support multiple pedigree relationships across sampling packs
Pairwise relatedness estimates ranged from 0 between putatively unrelated wolves to 0.484 between first-order relatives (such as parent-offspring and full sibling pairs). Although relatedness within packs was significantly higher than between packs (t test, t = 5.95, p < .001; Table S3)

| Different body sites exhibit distinct microbial communities
Alpha and beta diversity significantly differed between body sites using multiple diversity metrics. Examination of observed ASVs (Kruskal-Wallis test; H = 33.83, df = 6, p < .001) revealed higher species richness at heavily haired body sites and lower species richness at mucosal and fecal sites (Figure 3a). This pattern was largely driven by significantly lower diversity in mucosal nostril and perianal swabs when compared to haired dorsal flank and axilla swabs (Benjamini-Hochberg corrected q-values < 0.05; Table S5).
We also observed significant differences in species equitability between body sites using Pielou's evenness (H = 48.144, df = 6, p < .001), with high equitability observed in fecal and lip communities, moderate equitability at haired sites, and low equitability in the nostril (Figure 3b). Pairwise comparisons returned significant differences between perianal swabs and other body sites (except feces and lip), and nostril swabs and other body sites (except feces; q-values < 0.05; Table S6).
Principal coordinate analysis revealed sample clustering by body site, with differences in centroid position and dispersion evident ( Figure 3c; Figure S4). These patterns were supported by hierarchical clustering of body sites using Euclidean distances, where perianal and fecal samples formed one group, mucosal lip and nostril samples formed another, and haired axilla, dorsal flank, and ear canal samples formed the third ( Figure S5). This evidence suggests that microbial species richness, evenness, presence, and abundance all contribute to the body site-specific patterns of diversity observed in YNP wolf microbiomes.
Single-factor PERMANOVA of dorsal flank samples (Table S7) similarly returned sampling pack as the variable explaining the most variance within skin microbial communities, although this result was only statistically significant using unweighted UniFrac F I G U R E 2 Pedigree relationships for the 24 wolves included in this study (shaded) and their parents. Colors correspond to sampling pack, and dashed lines connect the same wolf in disparate parts of the pedigree distances (R 2 = .256, df = 3, p = .007). Field season was also significantly associated with species presence (unweighted UniFrac, R 2 = .202, df = 2, p = .001) in single-factor analyses. Multifactor PERMANOVA containing all six explanatory variables (Table 2;   Table S8) Table S9). Similar results were obtained with the 86,545 SNP dataset for perianal samples, and no statistically significant results were obtained for dorsal flank samples (Table S9).

| D ISCUSS I ON
In the present study, we characterized host-associated bacterial communities inhabiting haired and mucosal body sites of gray wolves sampled in situ at Yellowstone National Park. Consistent with our hypotheses, we reported distinct microbial communities colonizing each body site, with haired sites hosting high alpha diversity and gut samples harboring the largest proportion of anaerobic bacteria. We additionally contextualized skin and gut microbiota within genetic, demographic, and environmental variables to identify the primary factors underlying microbial community presence and abundance.
We found that social environment and host genetics were most consistently associated with microbial composition at both body sites.
We additionally reported evidence of body condition influencing gut microbiota, and coat color influencing skin microbiota. These results advance the field of microbial natural history and contribute to our understanding of the evolutionary history, ecology, and conservation of gray wolves and their associated microbes.
Examination of resident phyla enabled high-level contextualization of wolf-associated microbiomes within Canidae (Table 3).
In addition to high-level taxonomy, examination of microbial diversity and genera unique to each body site provides finer-scale insights into the form and function of commensal microbes. As seen in humans (Grice & Segre, 2011 , 2020), and may result from increased exposure to exogenous environmental factors. In contrast, lower species richness in mucosal sites suggests colonization by microbes specifically adapted to local conditions, as seen in sebaceous human skin (Grice & Segre, 2011).
At the genus and family levels, high proportions of aerobic bacteria in wolf nostrils and anaerobic bacteria in wolf guts suggest that oxygen exposure may contribute to differences between respiratory and gastrointestinal microbiota (Dickson & Huffnagle, 2015;Huffnagle et al., 2017;Jalili-Firoozinezhad et al., 2019;Thursby & Juge, 2017).
These microbes may further function in body site-specific processes. For example, dominant nostril families Moraxellaceae, Mycoplasmataceae, and Neisseriaceae have been isolated from the noses of domestic dogs and cats, with imbalances associated with nasal disease Tress et al., 2017). Similarly, dominant perianal genera exhibit numerous digestive and healthpromoting properties. For example, Bacteroides and Clostridium spp. function in bile acid metabolism and digestion of animal proteins and saturated fats (Deng & Swanson, 2015). Although associated with gastrointestinal disease in humans (Hussan et al., 2017), Fusobacterium spp. likely contribute to the breakdown of amino acids in healthy dogs and other mammalian carnivores (Pilla & Suchodolski, 2020;Vázquez-Baeza et al., 2016;Vital et al., 2015).

TA B L E 3
The different gray shades serve as a simplified heat map-the lower the rank (indicating higher abundance), the darker the shade of gray. function in lipid metabolism (Yang et al., 2020). The short-chain fatty acids produced by these genera (e.g., butyrate) further possess anti-inflammatory properties that may promote overall gut health (Tomova et al., 2019). Given significant differences between body sites, we independently identified factors underlying skin and gut microbiota, respectively. Sampling pack emerged as the variable most consistently associated with microbial presence and abundance at both body sites, explaining 31.7%-42.3% of the variance in gut samples and 20.8%-25.6% of the variance in skin samples. This result is consistent with social microbiome theory (Grieneisen et al., 2017;Sarkar et al., 2020) and gray wolf social ecology (Mech & Boitani, 2003;Stahler et al., 2020), as pack members cooperatively share territory, prey, social contacts, and microbes. Similar patterns have been ob- To disentangle the role of relatedness in these patterns, we leveraged biobanked samples and the highly resolved YNP wolf pedigree  to evaluate microbial similarity between first-order relatives, distant relatives, and unrelated dispersers with short (e.g., 1108M and 1107M) and long (e.g., 1047M) histories of dispersal from their natal packs. We found significantly positive correlations between host genetic and microbial distances at both body sites. However, after controlling for shared sampling pack, this result only remained significant for gut samples. This suggests that host genetics influence gut microbes more than skin microbes, perhaps due to differences in environmental exposure or maintenance of specialized microbiota (Grice & Segre, 2011;Thursby & Juge, 2017).
Similar patterns have been observed in humans, where cohabitation exerted a stronger effect on the skin microbiome compared to fecal microbiota (Song et al., 2013), as well as in great spotted cuckoos (Clamator glandarius) and magpies (Pica pica), where shared diet and environment were predictive of esophageal but not cloacal microbiota (Lee et al., 2020).
We observed further differences in the secondary factors underlying gut and skin microbial communities. For example, body condition significantly explained 30.4% of variance in gut microbial species abundance (N.B., this value reduced to 13.3% when controlling for sampling pack). This follows previous studies linking body condition, health, and disease to microbial composition in a variety of contexts, including clinical medicine (Gupta et al., 2020), veterinary medicine (Bradley et al., 2016;Rodrigues Hoffmann et al., 2014), and wildlife rehabilitation and conservation (DeCandia et al., 2019. We additionally found that coat color significantly explained 7.20% of variance in skin microbial species abundance when controlling for sampling pack. This presented a compelling preliminary result, as the mutation underlying melanism in gray wolves encodes for a beta defensin protein known to regulate microbiota and immune processes at the skin barrier (Anderson et al., 2009;Candille et al., 2007;Meade & O'Farrelly, 2019). Consequently, this mutation may lead to characteristic differences in the skin microbiota of black and gray wolves: an exciting frontier for further study.
Characterizing body site-specific microbiota and the factors underlying them is an important objective within evolution and ecology. While sampling wildlife populations in situ present logistical challenges, these analyses are critical for informing our understanding of the evolutionary history, ecology, and conservation of hosts and their associated microbes. Regarding evolution and ecology, we found that the core microbiota of YNP wolves mirrors numerous species within Canidae, likely due to shared phylogenetic history and characteristics such as carnivory. We additionally found that different body sites host distinct microbiota, with functions and products directly related to physiological processes (e.g., aiding digestion) and host health (e.g., regulating inflammation). We further found that social environment and host genetics influence skin and gut microbiota, with evidence supporting body condition and coat color as secondary factors influencing body site-specific microbiota.
Considered together, these results provide important baseline information for the long-term conservation of gray wolves and related species (Trevelline et al., 2019;West et al., 2019). We can now monitor in situ populations for the emergence of novel pathogens or microbial imbalances that may contribute to morbidity and mortality on the landscape (DeCandia et al., 2018). We can further adjust the diet (e.g., increase meat consumption) and social environment (e.g., enable group living) of captive-managed gray wolves to promote more natural microbiota. This study represents an important step toward these conservation goals and calls for similar studies to be conducted in additional host-microbe systems in captivity and the wild. Host characteristics such as evolutionary history, genetics, ecology, demography, and environment shape microbiota in innumerable, nuanced ways-and microbiota, in turn, shape their hosts. By characterizing the microbiota inhabiting diverse wildlife systems across body sites and environmental contexts, we can elucidate common patterns and processes that deepen our understanding of these relationships and contribute to successful wildlife conservation efforts.

ACK N OWLED G M ENTS
We would like to thank all of the Yellowstone Wolf Project technicians who collected observational data on YNP wolves. We would also like to thank the Smithsonian Conservation Biology Institute (SCBI) for providing salary support for ALD during data analysis and manuscript preparation, and the SCBI microbiome journal club for insightful conversations. This work was supported by Yellowstone