Fine-Scale Differentiation in Diet and Metabolomics of Small Mammals Across a Sharp Ecological Transition

Sharp ecological boundaries often present animals with abrupt transitions in various resources, including the availability of food, potentially creating strong selective gradients. Yet little is known concerning how animals respond to abrupt shifts in resources, especially when gene flow may limit local adaptation. Here, we used DNA metabarcoding and untargeted metabolomics of fecal samples to begin characterizing the foraging ecology of two closely related rodents (Neotoma bryanti and N. lepida) across a sharp ecological boundary. We find abrupt transitions in diet and metabolomic signatures that coincide with the rapid habitat transition. Further, we show that individuals have habitat-specific diets that are dominated by distinctly toxic plants that likely require different metabolic processing. Our approach allows detailed characterization of diet and metabolic processing of food plants, which can provide evolutionary ecologists and wildlife biologists much needed insight into the nutritional and physiological ecology of the systems they study and manage.


INTRODUCTION
Strong selection across sharp ecological boundaries can produce or maintain patterns of adaptive divergence over small spatial scales, even in the presence of gene flow (Antonovics, 2006;Grahame et al., 2006;Rosenblum, 2006, reviewed in Richardson et al., 2014. Sharp ecotones provide spatially proximate but selectively distinct sets of environmental factors to which organisms may respond, thereby potentially initiating species divergence or influencing the evolutionary trajectory of species coming back into secondary contact after having diverged in isolation (Smith et al., 1997;Schluter, 2001). For animals, one of the main environmental factors to which they must respond at sharp habitat boundaries or in heterogeneous environments is food availability. Indeed, some of the best examples of ecological adaptation are the result of diet-related selection, including bill variation in Darwin's finches (Grant and Grant, 2006), gill raker morphology of stickleback fishes (Schluter and McPhail, 1992), and jaw morphology in salamanders (Adams and Rohlf, 2000). For some herbivores, food consumption is ultimately a chemically-mediated ecological interaction wherein animals must balance the acquisition of nutrients while minimizing exposure to potentially toxic plant secondary compounds . Some of the primary diet-related adaptations in these herbivores, then, are likely associated with how these animals protect themselves metabolically from overexposure to plant toxins (Skopec et al., 2008;Forbey et al., 2018). Identifying the underlying phenotypic plasticity and evolutionary potential (Sorensen et al., 2005;Magnanou et al., 2009;Gompert et al., 2015;Oh et al., 2019) of the traits that allow metabolic processing of plant toxins is central to understanding dietrelated ecological adaptation in herbivores. Further, identifying how animals respond to changes in exposure to plant secondary compounds, spatially and temporally, will be critical in predicting how animals will respond to climate change, and the associated shifts in community composition and food availability (Estrada et al., 2016;Beale et al., 2018;Forbey et al., 2018).
Here, we take advantage of the overall ecological transition between the relatively mesic southeastern slopes of the Sierra Nevada Mountains and the arid western valleys of the Mojave Desert (Figure 1) to understand how small mammal herbivores of the genus Neotoma (woodrats) respond to fine-scale shifts in plant availability and toxin exposure. At the western edge of the Kelso Valley of California there is a sharp habitat transition (within tens of meters) between the mesic woodland uplands, or "hill" habitat, and the starkly arid desert basin, or "flats" habitat immediately adjacent to the base of the hill (Figure 1; Shurtliff et al., 2014). At this site, the relatively mesic-adapted Neotoma bryanti (Bryant's woodrat) primarily occupies the boulder-strewn hill habitat, while the relatively aridadapted Neotoma lepida (Desert woodrat) primarily occupies the flats habitat Shurtliff et al., 2014). Although the adjacent hill and flats habitats are well within the dispersal capabilities of individual woodrats, multi-year mark-recapture studies have shown that there is limited short or long-term movement between the two habitats (Shurtliff et al., 2014). Despite the overall spatial segregation of the two species, approximately 13% of the woodrats at the site have hybrid ancestry, ranging from first generation hybrids to highly backcrossed individuals in both parental directions that are largely co-distributed with their predominant parental genomic type (i.e., backcrossed-lepida are primarily found in the flats and backcrossed-bryanti are found almost solely on the hill; Shurtliff et al., 2014).
We hypothesize that the strong spatial segregation of the species and their respective backcrosses at this site is at least in part due to differential adaptations in metabolic processing of available plants in these distinct, yet adjacent, habitats. As such, we predicted that woodrats in the flats habitat would consume a diet distinct from woodrats on the hill. Further, many of the plants specific to either habitat are phytochemically distinct and well-defended including sagebrush (Artemesia tridentata), rabbitbrush (Ericameria sp.), and desert almond (Prunus fasciculata) in the flats; and California coffeeberry (Frangula californica) and pine (Pinus sp.) on the hill. As such, we sought to develop chemical profiling methods to identify fecal metabolites that may be associated with dietspecific metabolic processes. Thus, in addition to predicting we would find habitat-specific diets, we also predicted that woodrats would exhibit habitat-specific metabolomic profiles as a result of their encounters with unique phytochemical toxins in their diets. This system provides a unique opportunity to investigate the potential for fine-scale differential plant consumption and possible distinct strategies for metabolic processing of plant toxins even in the face of hybridization and gene flow between distinct habitats.

Field Sampling
Adult woodrats build and defend a large stick nest that they occupy solely, and each nest has an obvious latrine area where fecal pellets accumulate. To begin identifying which plants woodrats are consuming in these habitats, in July 2013, we collected fecal pellet clusters from 28 unique woodrat houses known to have been previously occupied by N. bryanti on the hill (n = 15) and N. lepida in the flats habitat (n = 13). We expected that this initial snapshot of the diet in these habitats would reflect among the most restrictive dietary conditions these animals face. That is, in July of any year at this desert site, plant availability is likely to be primarily restricted to perennial shrubs and trees, but because 2013 was a particularly severe drought year, we expect that few annual forbs were available in the weeks leading up to our sampling. We collected samples at nests within a 150 m radius of where the base of the hill meets the flats (Figure 1) in order to capture the abrupt habitat transition. Beyond being associated with obvious woodrat nests, all fecal pellets collected were the characteristic shape and size of adult Neotoma (Smith et al., 1995). Although samples were clearly associated with a single nest and our genotyping did not show evidence of multiple individuals within a sample (see below), we cannot be completely certain that all pellets taken from a nest were left by the same individual, and thus, these should be considered community samples representative of the surrounding habitat. We stored fecal samples in paper envelopes and allowed them to continue to dry at ambient temperature until processed. All plant metabarcoding and chemical analyses were conducted within 1 year of sample collection.

Woodrat Genotyping
We classified individual samples to species by extracting DNA from three fecal pellets per house using the QIAamp DNA Stool kit (Qiagen, Inc.). We followed manufacturer instructions with the following modifications: (1) we allowed fecal pellets to soak in the first buffer for 4 h, then pipetted out the buffer to start the first step of the protocol, and (2) the final elution step was done with 100 µl of molecular grade water, and this elute was passed through the filter twice to increase DNA FIGURE 1 | Diet composition of woodrats occupying the "flats" and "hill" habitats of a sharp ecological transition that exists at the Whitney Well site of Kelso Valley, California. The flats habitat is predominantly occupied by N. lepida while the hill is largely occupied by N. bryanti. Diet is represented as the relative read abundance of plant DNA found in fecal samples collected at woodrat nests in each habitat (N hill = 13, N flats = 12). A list of all plants, the range of relative read abundance values, and frequency of occurrence in each habitat are provided in Supplementary  concentration. We genotyped samples at a suite of microsatellite loci developed by Sousa et al. (2007) and following lab methods as described in Coyner et al. (2015). Species identity was made by combining the samples of this study with a larger dataset of known N. bryanti and N. lepida from the study site and running STRUCTURE (Pritchard et al., 2000) at K = 2 as in Shurtliff et al. (2014).

Fecal DNA Barcoding and Diet Composition
Fecal pellets were sent to SPYGEN (Inc.) where plant DNA was extracted from 8 to 10 individual pellets (∼10 mg) per pellet cluster using the DNeasy Mini Stool Kit (Qiagen, Inc.) following manufacturer instructions. A segment of the trnL gene was amplified using the universal primers gh that typically result in fragments between 10 and 143 bp depending on the plant species (Taberlet et al., 2007;Pompanon et al., 2012). DNA amplification was repeated 4 times per sample and mock extractions and amplifications were carried through each experiment to monitor for contamination. Amplified products were titrated by capillary electrophoresis (QIAxcel; Qiagen GmbH, Hilden, Germany) and purified using the MinElute PCR purification kit (Qiagen GmbH, Hilden, Germany), then titrated again using capillary electrophoresis. Purified PCR products were pooled in equal volumes. Library preparation and sequencing were performed at Fasteris (Geneva, Switzerland). Libraries were prepared using the TruSeq NanoDNA genomic kit (Illumina, San Diego, CA, United States) and paired-end sequenced (2 × 100 bp) on an Illumina MiSeq (Illumina, San Diego, CA, United States) using the Paired-end MiSeqReagentKit V2 (Illumina, San Diego, CA, United States) following manufacturer's instructions. Raw sequences were demultiplexed and quality filtered and eliminated from the dataset if they had less than 98% identity to plant sequences held in GenBank (following Valentini et al., 2009). Putative OTUs (operational taxonomic units) were established based on clustering sequences with > 98% sequence similarity, and the most commonly recovered sequence was designated as representative for that OTU; this led to 33 OTUs represented by at least 50 reads across the entire dataset. Plant identification from these representative reads was made by conducting a BLAST search of the NCBI standard databases 1 using a blastn search of "somewhat similar sequences" and inspecting all taxa with > 98% sequence identity. When only a single species had 100% sequence identity and was known to occur at the site from our own collections or from the region based on the CalFlora database 2 , we identified the read as representing that plant species. When multiple species of the same genus had 100% identity with the searched read, we narrowed the blastn list to those species known from the region based on the CalFlora database or from our own collections at the site, and only report the genus level since multiple species could be represented by the sequence we recovered. Finally, when multiple species within multiple genera known from the region or our own collections matched our read at 100%, we report the family level (Supplementary Table 1). To estimate diet composition we established the relative read abundance (RRA) of each plant within each woodrat sample, and then averaged across all hill or all flats samples to estimate habitat-based diet composition. Here, we assume that RRA is proportional to the amount of a food item consumed (Bergmann et al., 2015;Kartzinel et al., 2015). While we recognize that RRA may only be a coarse index of quantity of food consumed, even with moderate recovery biases, RRA has been shown to provide an accurate view of dietary composition at the population level (Deagle et al., 2019). We also report the frequency of occurrence (FOO); a measure that only relies on the presence/absence of plants rather than their abundance (Deagle et al., 2019). Finally, while we have confirmed that these primers amplify the most common shrub species at the study site, we have yet to conduct thorough digestion studies (e.g., Willerslev et al., 2014), so our results should be viewed as a coarse, overall view of diet composition across this study site at the time of sampling.

Fecal Metabolite Analyses
Operating with no prior knowledge of what woodrats consume at this site, we applied an untargeted metabolomics approach to generate chemical profiles from the same fecal samples used in the DNA diet analysis (see above). Our approach is designed to generate chemical profiles from small samples (fecal pellets) in order to quantify chemical variation at the individual level, allowing for direct comparison between individual plant consumption and chemical profile. We optimized a protocol adapted from Santos Pimenta et al. (2014), which entailed extracting 200 mg dry weight from each of 28 fecal samples in 2 mL of a 50:50 KH 2 PO4 buffer:CD 3 OD. We analyzed extracts using both proton nuclear magnetic resonance ( 1 H-NMR) spectroscopy and high-pressure liquid chromatography coupled with mass spectrometry (HPLC-MS) to achieve untargeted metabolite profiles. HPLC-MS analysis is more sensitive than 1 H-NMR and can detect compounds that occur at very low concentrations, but is more dependent on the specific ionization efficiencies of the molecules under analysis when compared to 1 H-NMR. To minimize the effect of the KH 2 PO 4 buffer on the HPLC-MS chromatograms, we added 3 mL H 2 O + 0.1% formic acid and lyophilized the aqueous suspension prior to obtaining extract masses. We then dissolved samples in methanol to approximately 1 mg/mL for HPLC-MS analysis. Samples were analyzed using an Agilent 1200 HPLC coupled to an Agilent G6230A high-resolution high-mass accuracy time-offlight (TOF) mass spectrometer with an electrospray ionization (ESI) source using the following settings: gas temperature: 325 • C; gas flow: 8 L/min; nebulizer: 35 psig; capillary voltage: 3500 V; fragmentor: 175 V; skimmer1: 65 V; octopole: 750 V; mass range: 100-3000 m/z. Chromatographic separation was accomplished using a Phenomenex Kinetex-EVO C18 2.6 µm column (2.1 × 100 mm, 100 Å) maintained at 40 • C. The mobile phase consisted of (A) H 2 O + 10 mM NH 4 OAc and (B) MeCN with an elution gradient of: 1-14 min, 5-100% B; 14-16 min, maintain 100% B; 16-17 min, 100-5% B; 17-20 min, maintain 5% B for column re-equilibration. The flow rate was 0.5 mL/min, and the sample injection volume was 2 µL. Finally, we analyzed extracts from several individuals of the dominant diet plants in each habitat. Approximately 200 mg dry weight from each individual was extracted in 3 mL of 50:50 CH 3 OH:H 2 O. Samples were vortexed for 30 s, sonicated for 15 min, and centrifuged for 15 min before filtering. An aliquot of each sample was diluted with H 2 O for HPLC-MS analysis as described above. The remaining liquid extract was dried down under nitrogen and reduced pressure. Samples were then reconstituted in 50:50 CD 3 OD:D 2 O for 1 H-NMR analysis.
We processed 1 H-NMR data using MestReNova software (Mestrelab Research, Spain). We aligned spectra from crude extracts using the solvent peak, then baseline and phase corrected, and binned every 0.04 ppm from 0.5 to 12 ppm using the global spectral deconvolution method. We removed solvent and water peaks prior to normalizing bins to a total area of 100. We ran a Mantel test using Spearman's rank correlation coefficients in R (mantel function in vegan, Oksanen et al., 2018;R Core Team, 2019) to compare the Bray-Curtis dissimilarity matrix of the chemical profiles of fecal samples to the dissimilarity matrix of the plants present in each fecal sample based on DNA data. We used Partial Least Squares Discriminant Analysis (plsda function in mixOmics, Lê Cao et al., 2016) to visualize the 1 H-NMR data. The HPLC-MS data were processed using Proteowizard (Chambers et al., 2012) and MZmine (Katajamaa et al., 2006) to obtain an aligned peak list. We used a network approach implemented in Cytoscape (Shannon et al., 2003) to visualize relationships between diet, habitat and fecal chemical profiles. We tested for differences in the chemical profiles of fecal samples from hill and flats habitats using a permutational multivariate analysis of variance of the Bray-Curtis dissimilarity matrices (function adonis in vegan package). For the HPLC-MS data, we also used Non-Metric Multidimensional Scaling (NMDS) for dimension reduction and data visualization. We tested whether the proportion of the dominant diet plants found through DNA analysis predicted chemical dissimilarity by projecting vectors on the NMDS plot that maximized the correlation between chemistry and proportion of the dominant diet plants (function envfit in vegan package). Finally, we validated these results using a linear regression model in which the proportion of the dominant diet plants in the sample predicted the first coordinate scores.

RESULTS
Using microsatellite genotyping, we confirmed that 12 of 15 hill samples were left by N. bryanti, and 10 of 13 flats samples were left by N. lepida. Six woodrat fecal samples did not return an adequate number of alleles to definitively assign them to one species or the other. Because the results we report below are qualitatively the same with and without the woodrat samples that we could not attribute to one species or the other, we report data from all 25 samples that yielded trnL data that passed quality filtering.
From plant trnL sequencing, we recovered 2,345,967 reads from 25 sets of fecal samples resulting in a final sample size of N hill = 13 and N flats = 12 because we eliminated two of the original 15 hill samples and one of the original 13 flats samples that returned fewer than 500 reads each. The trnL dataset provided evidence of 33 unique taxa that included resolution to the species (19), genus (7), or family (7) level (Supplementary Table 1). All taxa recovered are well-known from the region (according to the CalFlora database) and most have been documented at the site and collected by our group (Q. Shurtliff and D. Nielsen, pers. comm.). Eleven plant taxa were recovered in diets from both the hill and the flats; another 11 taxa were only recovered from the flats; and the final 11 taxa were only recovered from the hill (Supplementary Table 1).
On the hill, 10 of 13 samples contained F. californica (frequency of occurrence (FOO) = 77%; Supplementary Table 1), while in the flats, 10 of 12 samples contained P. fasciculata (FOO = 83%). In contrast, one of 13 hill samples showed evidence of P. fasciculata (FOO = 7.7%), while one of 12 flats samples showed evidence of F. californica (FOO = 8.3%). The only other plant with a frequency of occurrence greater than 50% in either habitat was E. umbellatum. This species of buckwheat was found in 8 of 13 samples from the hill habitat (FOO = 62%). Only one sample from the flats habitat showed consumption of buckwheat (FOO = 8.3%), and this was the same sample that showed consumption of F. californica in the flats.
We were able to obtain HPLC-MS data from N hill = 10 and N flats = 10 fecal samples for which we present DNA data, but because 1 H-NMR requires higher sample concentration than HPLC-MS, we were only able to obtain 1 H-NMR data from N hill = 7 and N flats = 8 samples. Using 1 H-NMR data, we found fecal samples have a significant correlation between their chemical profiles and plant RRA (Mantel statistic r = 0.19, p = 0.04). PLS-DA of 1 H-NMR data showed that hill and flats animals had distinct metabolomic profiles (Figure 2A, R2 = 0.88, Q2 = 0.63), although further analyses are needed to interpret the specific 1 H-NMR signal loadings of the PLS-DA. HPLC-MS data showed that there were many metabolites that were ubiquitous across all samples, but also revealed several habitat-specific metabolites [ Figure 3, habitat effect on chemistry F (1 , 19) = 2.61, p < 0.01]. Consistent with the 1 H-NMR data, we found a significant correlation between chemistry and plant DNA found in fecal samples (Mantel statistic r = 0.26, p < 0.01). We found that the proportion of P. fasciculata and F. californica in the fecal samples was significantly correlated with chemical similarity (r 2 = 0.34, p < 0.05 and r 2 = 0.38, p < 0.05, respectively) in opposing directions along the first coordinate (NMDS1) of the NMDS plot. In a linear regression, we found that the proportion of F. californica in the diet significantly predicted chemistry [R 2 = 0.36, F (1 , 18) = 9.92, p < 0.01] and was positively related to NMDS1 with a slope coefficient of 0.58. In comparison, the proportion of P. fasciculata in the diet significantly predicted chemistry [R 2 = 0.33, F (1 , 18) = 8.80, p < 0.01], but was negatively related to NMDS1 with a slope coefficient of −0.58. Finally, we show that F. californica and P. fasciculata are chemically diverse plants that have distinct chemical profiles (Figure 4). Further analyses are required for full chemical characterization of these plants, including across seasons, individuals and tissue types. Nonetheless, our data show that F. californica from our study site contains a clear peak with a characteristic mass associated with emodin, a common anthraquinone (Figure 4 and Supplementary Figure 1). In P. fasciculata from our site we found evidence of masses and fragmentation that are characteristic of prunasin, a cyanogenic glycoside (Figure 4 and Supplementary Figure 1).

DISCUSSION
Our study demonstrates fine-scale diet and metabolomic differentiation across a sharp ecological boundary. Using genetic analysis of fecal samples we show that woodrats occupying each side of this sharp ecotone largely consumed distinct sets of plants during our sampling period. Frangula californica was the most common diet item in the hill habitat (FOO = 77%, RRA = 52%) while Prunus fasciculata was dominant in diets of woodrats in the flats habitat (FOO = 83%, RRA = 59%). Additional sampling across seasons and years will be needed to identify the importance of Frangula and Prunus in the dietary ecology of these animals. Nonetheless, that woodrats can consume large amounts of either of these plants for any duration of time is surprising because both are known to contain compounds that are potentially highly toxic to herbivores. From the literature we know that P. fasciculata contains cyanogenic glycosides that are transformed into hydrogen cyanide upon mastication and later in the digestive process (Santamour, 1998), while F. californica contains anthraquinones, which cause damage to the intestinal lining resulting in a strong laxative effect (Surh et al., 2013). We found evidence that the plants at our study site had characteristic HPLC-MS fragments for both cyanogenic glycosides in P. fasciculata (i.e., prunasin) and anthraquinones in F. californica (i.e., emodin). Holistically, HPLC-MS data of these plants reflect distinct chemical compositions that likely require different detoxification strategies (see below). Therefore, we anticipated discovering distinct metabolic products in the feces of the woodrats consuming these plants.
Our chemical analyses did, indeed, reveal distinct chemical signatures in the feces of animals consuming these different diets. Specifically, woodrats on the hill, that have diets containing a large amount of F. californica, had several unique fecal metabolites not shared with woodrats from the flats that consume diets dominated by P. fasciculata, and vice versa (Figure 3). Moreover, while both hill and flats animals eat a variety of plants, the proportion of P. fasciculata and F. californica consumed based on DNA evidence is the primary driver of differentiation between the fecal metabolomic profiles of these animals ( Figure 2B). The habitat-specific chemical profiles we document here are a composite from various sources including unaltered plant chemistry (undigested plant material) as well as plant material chemically modified by woodrat and microbial digestive and detoxification processes. Further analyses, including controlled feeding trials with individual plants, will be needed to characterize the unique metabolites we have discovered here. Likewise, further sampling of woodrat tissues/products including blood and urine will be needed to capture a more complete array of dietary compounds from these plants and their post-ingestion metabolites (Dunn et al., 2011;Forbey et al., 2018). Nonetheless, these fecal metabolomic profiles reflect the distinct chemistries of the habitat-specific diets in this system, along with how hill and flats woodrats are metabolically processing these diets.
Other woodrat populations are known to consume large quantities of highly toxic plants including creosote bush (Larrea tridentata; Magnanou et al., 2009), juniper (Sorensen et al., 2005), oak leaves (Skopec et al., 2008), and mesquite (Smith et al., 2014). Our data suggest that the woodrats occupying the hill habitat (predominantly N. bryanti) have developed or evolved the capacity to consume large quantities of F. californica, while woodrats on the flats (predominantly N. lepida) have developed or evolved the capacity to consume large amounts of P. fasciculata. Our study adds to the growing list of chemically well-defended plants that members of the genus Neotoma appear FIGURE 3 | Tripartite network relating plant diet (Frangula californica and Prunus fasciculata) to individual woodrat fecal samples from the hill (H) and the flats (F) habitats and the metabolites (m/z features) identified through LC-MS of their respective fecal samples. The metabolites cluster into three groups, those that are shared between all rats (gray), those unique to hill rats consuming F. californica (green) and those unique to flat rats consuming P. fasciculata (pink).
capable of consuming in large quantities, and in some cases, on which they locally specialize.
Detoxification of F. californica and P. fasciculata has not been studied in woodrats, but data from laboratory rats and mice allow us to make initial predictions of metabolic processing of the compounds held in these plants. Anthraquinones produced by F. californica have a strong laxative effect and, in other mammals, detoxification of these compounds is at least partly mediated in the liver by cytochrome P450 1A2 (Mueller et al., 1998). In contrast, the hydrogen cyanide that animals are exposed to when consuming P. fasciculata is partly detoxified by the enzyme rhodanese, which is heavily expressed in the liver (Dooley et al., 1995). It is likely that woodrats at least partly use these mechanisms to protect themselves from toxin exposure when consuming F. californica and P. fasciculata, respectively. In addition to their own detoxification pathways, the microbiome may also be involved in helping woodrats detoxify these compounds. In the case of creosote, the microbiome facilitates the woodrat's capacity to detoxify this plant (Kohl et al., 2014) suggesting that the microbiome may play a role in the detoxification of other plants woodrats consume (reviewed in Kohl and Dearing, 2016). Finally, the degree to which woodrats are exposed to either prunasin or emodin will also be determined by the interaction of these compounds with others that occur in these plants (Figure 4 and Supplementary Table 2) and those in other plants they consume simultaneously. It is intriguing to note that 12 of 14 individuals (86%) that consumed California coffeeberry (F. californica) also consumed sulpher buckwheat (E. umbellatum). This apparent co-occurrence could suggest that the combined chemistry of these plants might synergistically/antagonistically provide added nutrition and/or protection from toxins (Caesar and Cech, 2019).
This system presents a case of particularly fine-scale differentiation in feeding ecology despite gene flow across these species and habitat boundaries. Although our snapshot of diets at this site may represent among the most restrictive dietary conditions these populations face, the strong selection associated with these periods may have profound effects on population persistence and the evolutionary trajectory of these populations. Additional experimental work is needed to determine if N. lepida and N. bryanti have evolved unique genetically-based mechanisms to metabolize the dominant compounds in either P. fasciculata or F. californica, or whether they deploy existing metabolic pathways. Further, it will be critical to assess whether these species and their hybrids are capable of metabolizing the plant they are not adapted or acclimated to consuming in their native habitat (i.e., hill animals eating P. fasciculata and flats animals eating F. californica). Because neither plant is very common in either habitat (although P. fasciculata is predominantly found in the flats and F. californica is largely restricted to the hill), the distribution of these plants could influence the spatial distribution of these animals and the opportunity for interspecific contact and mating (prezygotic isolation, e.g., Caillaud and Via, 2000). Further, if only certain interspecific genomic or genomic-microbial combinations can successfully metabolize these diets, this may influence patterns of hybrid survival (postzygotic isolation) and genomic introgression across this system. Importantly, because there is interspecific gene flow across the habitat boundary, there is the potential for introgression of metabolically advantageous alleles and their integration into largely heterospecific genomes. Given this, the microsatellite species identifications presented here should be viewed as a coarse view of genome composition, because these animals could all have some heterospecific alleles in their genomes. Overall, systems such as this have the potential to provide novel insights into how a spectrum of mammalian genomes across a hybrid zone interact with microbial genomes to face the ecological challenges posed by their environment, including the plants they consume.

FUTURE DIRECTIONS
One of the most pressing challenges vertebrate populations will face as they respond to a changing climate is food acquisition. As species respond to changing conditions, shifts in plant phenology and distribution (Parmesan and Hanley, 2015) may reduce the availability of preferred food items for herbivores. In addition, environmental change (e.g., elevated CO 2 and temperature) is expected to dramatically modify plant chemistry, often increasing potentially toxic plant secondary compounds (Zvereva and Kozlov, 2006;Jamieson et al., 2017). Further, as temperatures increase, the rates at which herbivores can metabolically process plant secondary compounds are expected to decrease (Kurnath et al., 2016;Beale et al., 2018), further exacerbating increases in plant toxicity (Forbey et al., 2013). There will be growing need for ecologists and wildlife managers to monitor the metabolic capabilities of wildlife populations as plants and animals respond to changing conditions.
The integrated approach we demonstrate here provides a non-invasive method to quickly gain insight into wild herbivore diet composition, the chemical profile of plants they are predominantly being exposed to, and the products of metabolism that indicate how they are metabolizing particular plants and their secondary compounds. The case study we present here could be classified as a "rapid assessment" to create an initial baseline of diet composition and metabolic processing. Thorough reconstruction of diet would require sampling across seasons and years, and ideally, repeated sampling of known individuals. If fecal samples can be collected directly from live-trapped individuals and immediately frozen, it is possible to recover the fecal microbiome (see Kohl et al., 2015) in addition to plant DNA and fecal chemistry, providing an integrated view of concerted change in diet, microbial community composition, and overall metabolic processing. In addition to sampling animal diet, plant chemistry should also be monitored through time at a site, although the method we describe here is intended to recover relatively stable chemical compounds and metabolites, and not volatile compounds. Overall, in this study system and others, additional work is needed to (1) refine primer sets and DNA databases so that fecal DNA read counts more accurately reflect herbivore diets, (2) identify specific metabolites and the metabolic processes they reflect, and (3) integrate this approach into wildlife monitoring plans to link diet-related physiology to population health and demography (Forbey et al., 2018). The application of these approaches will require close collaboration between ecologists and wildlife managers, and will result in unprecedented insight into how these and other herbivores interact with their environment at a highly mechanistic level.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI GenBank (accession MT640105-MT640142).