Differences in dietary composition and preference maintained despite gene flow across a woodrat hybrid zone

Abstract Ecotones, characterized by adjacent yet distinct biotic communities, provide natural laboratories in which to investigate how environmental selection influences the ecology and evolution of organisms. For wild herbivores, differential plant availability across sharp ecotones may be an important source of dietary‐based selection. We studied small herbivore diet composition across a sharp ecotone where two species of woodrat, Neotoma bryanti and N. lepida, come into secondary contact with one another and hybridize. We quantified woodrat dietary preference through trnL metabarcoding of field‐collected fecal pellets and experimental choice trials. Despite gene flow, parental N. bryanti and N. lepida maintain distinct diets across this fine spatial scale, and across temporal scales that span both wet and dry conditions. Neotoma bryanti maintained a more diverse diet, with Frangula californica (California coffeeberry) making up a large portion of its diet. Neotoma lepida maintains a less diverse diet, with Prunus fasciculata (desert almond) comprising more than half of its diet. Both F. californica and P. fasciculata are known to produce potentially toxic plant secondary compounds (PSCs), which should deter herbivory, yet these plants have relatively high nutritional value as measured by crude protein content. Neotoma bryanti and N. lepida consumed F. californica and P. fasciculata, respectively, in greater abundance than these plants are available on the landscape—indicating dietary selection. Finally, experimental preference trials revealed that N. bryanti exhibited a preference for F. californica, while N. lepida exhibited a relatively stronger preference for P. fasciculata. We find that N. bryanti exhibit a generalist herbivore strategy relative to N. lepida, which exhibit a more specialized feeding strategy in this study system. Our results suggest that woodrats respond to fine‐scale environmental differences in plant availability that may require different metabolic strategies in order to balance nutrient acquisition while minimizing exposure to potentially toxic PSCs.

4. Neotoma bryanti and N. lepida consumed F. californica and P. fasciculata, respectively, in greater abundance than these plants are available on the landscapeindicating dietary selection. Finally, experimental preference trials revealed that N. bryanti exhibited a preference for F. californica, while N. lepida exhibited a relatively stronger preference for P. fasciculata. We find that N. bryanti exhibit a generalist herbivore strategy relative to N. lepida, which exhibit a more specialized feeding strategy in this study system. 5. Our results suggest that woodrats respond to fine-scale environmental differences in plant availability that may require different metabolic strategies in order to balance nutrient acquisition while minimizing exposure to potentially toxic PSCs.

| INTRODUC TI ON
Ecotones are characterized by spatial transition in environmental variables that can create selective gradients that generate or maintain diversity (Smith et al., 1997). When sharp abiotic gradients support the establishment of spatially proximate but distinct vegetation communities (Walker et al., 2003), animals must respond to abrupt spatial transitions in abiotic and biotic resources. Such spatially proximate, yet dissimilar selective environments have the potential to generate or reveal the ecological adaptations or forms of phenotypic plasticity that permit species to exist in disparate environments (Ghalambor et al., 2007;West-Eberhard, 1989).
At sharp environmental transitions, one of the primary challenges facing herbivores is the abrupt transition in food plant availability. For herbivores, space use and movement across ecotones is largely determined by the distribution of plants that allow acquisition of adequate nutrition, while minimizing exposure to toxic plant secondary compounds (PSCs; Dearing et al., 2000Dearing et al., , 2005Freeland & Janzen, 1974;Westoby, 1978). Mammalian herbivores have evolved numerous behavioral and physiological adaptations to maximize nutrition while minimizing toxin exposure including regulation of liver detoxification enzymes (Malenke et al., 2012), decrease in metabolic rate and physical activity when exposed to dietary PSCs (Sorensen et al., 2005b), and maintenance of a microbiome that facilitates nutrient acquisition and detoxification (Kohl et al., 2014). Mammalian herbivores may also diversify their diets to minimize overexposure to, or neutralize, toxins present (Iason & Villalba, 2006). Based on the degree to which mammalian herbivores modify their diets either spatially or temporally, they can be classified along a continuum of foraging strategies from generalist to specialist consumers (Shipley et al., 2009).
When mammals consume toxic plants they are not adapted to, they suffer energetic consequences that can lead to rapid weight loss and lowered body condition (Mangione et al., 2004;Sorensen et al., 2005aSorensen et al., , 2005b. Given these consequences, we would expect mammalian herbivores to develop dietary preferences for plants with which they are familiar and which they can efficiently digest (Partridge, 1981). Hence, for herbivorous mammals, distinct vegetation communities across sharp ecotones may produce spatial variation in selection that leads to or reinforces distinct dietary preferences, which may in turn determine fine-scale space use and a range of intra and interspecific interactions (Nosil et al., 2005;Via, 1999;Via et al., 2000).  Patton et al., 2007), and while they are largely spatially segregated between the two adjacent habitats, they do occasionally hybridize. These hybridization events lead to approximately 14% of individuals across the study site having hybrid ancestry, with backcrossing and introgression in both parental directions (Shurtliff et al., 2014;J. P. Jahner, T. L. Parchman, M. D. Matocq, unpublished data). Previous diet analyses (Matocq et al., 2020;Shurtliff et al., 2014) suggest that N. bryanti and N. lepida consume distinct diets in the hill and flats, respectively. As such, this system offers an opportunity to investigate dietary choices across a sharp ecotone, as well as the potential role of dietary differences in limiting interspecific contact and hybridization.
Here, we sought to further characterize the degree to which dietary composition and preference differ between pure N. bryanti and N. lepida in their respective native habitats, and to uncover the potential ecological correlates maintaining species differences in diet across this ecotone. We integrate both field and laboratory studies to ask the following questions: (1) Do N. bryanti and N. lepida maintain distinct diets across this sharp ecotone in both wet and dry seasons, and in wet and dry years? (2) Do these species consume plants in the wild in proportion to their availability in the habitat, or do they exhibit selection/preference for particular plants? (3) When given a choice in experimental trials, do woodrats exhibit the same dietary preferences as exhibited in field-collected samples? (4) To what degree are plant nutritional content and plant secondary compounds correlated with dietary preferences? To address these questions, we quantify diet preferences in the wild using high-throughput sequencing of the chloroplast trnL intron from woodrat fecal samples collected across the ecotone. We further examine these apparent patterns of preference by conducting an experimental choice trial.
To understand the underlying drivers of fine-scale diet differentiation in this system, we place these dietary preferences within the context of availability of these plants on the landscape, the plant secondary compound composition of these plants, and their nutritional quality. Our study provides a well-developed example of fine-scale diet differentiation in mammalian herbivores-differences across an ecotone that are maintained between the species in both wet and dry conditions.

| Study system
The study site is located in Kelso Valley, Kern Co., California, where N. bryanti and N. lepida meet and hybridize at the southern end of the Sierra Nevada mountain range (35°25′45N, K E Y W O R D S adaptation, detoxification, herbivore, hybridization, Neotoma, toxin tolerance, woodrat 118°15′2W). The mesic "hill" habitat sharply transitions to the xeric "flats" habitat ( Figure 1), and both parental species and hybrids can be found across a span of as little as tens of meters. The total area of the study site is approximately 50 hectares, approximately centered at the base of the hill (Figure 1). We conducted vegetation surveys in 27 plots (hill = 16, flats = 11) to estimate the abundance of the most common shrubs and trees (details in Supporting Information).

| Woodrat species identity
We identified individuals as N. bryanti or N. lepida using microsatellite loci previously developed for these species (Sousa et al., 2007).
For animals included in the preference trials (see below), we obtained ear biopsies from each individual and conducted DNA extraction, amplification and scoring of microsatellite loci as described in Coyner et al. (2015). We established species identity by conducting a Bayesian assignment test as implemented in STRUCTURE (Pritchard et al., 2000;Falush et al., 2003) at K = 2 as in Shurtliff et al. (2014) and used q lepida values >90% to assign individuals to N. lepida and q lepida values <10% to assign individuals to N. bryanti. To confirm the species identity of individuals included in the fecal metabarcoding, we used the same genotyping approach, but started with the gDNA extractions used for trnL sequencing (see below) and performed three replicate PCRs per sample. after 3-4 nights, we collected fresh, adult-sized fecal pellets. It is important to note that woodrat houses are solely occupied by one adult woodrat, and these animals are highly territorial, so there is limited chance that more than one woodrat contributed to the fresh fecal pellets we collected. We placed pellets into coin envelopes to F I G U R E 1 Panel (a) depicts the study site where the mesic hill transitions to the xeric flats. Photo taken from the north looking south. Black star in inset map represents approximate location of the study in Kelso Valley, California. Panels (b) and (c) depict habitat of the flats and hill habitats, respectively. Inset photo of woodrat is Neotoma lepida

| Fecal metabarcoding
dry, and stored them long-term at −20°C. We submitted samples to

| Diet composition
We used read counts of all identified plants to calculate Shannon diversity for diets of N. bryanti and N. lepida and performed a twosample t-test in R to compare diversity in diet composition (R Core Team, 2016). We used read counts to determine if diets between the two species were distinct by performing a PERMANOVA using Bray-Curtis distances with the adonis function in the vegan package (Oksanen et al., 2013;R Core Team, 2016).
To estimate individual and population-level (i.e., species) consumption of particular plants, we used both frequency of occurrence In order to take individual variation into account in estimates of population-level consumption, we used a hierarchical Bayesian approach implemented in R using the bayespref package (Fordyce et al., 2011) to estimate population-level consumption of the 5 most common plants identified in woodrat diets, which comprised ~80%-90% of total reads (Tables 1 and S1-S3). We pooled the remaining read counts from all other plant taxa into an "other" group. Rather than relying simply on RRA (as described above) to infer relative degree to which plants are consumed, this hierarchical Bayesian approach incorporates individual variation in our population-level consumption estimates (Fordyce et al., 2011;Forister et al., 2013).
We used raw read count data to run models. Raw read counts were not normally distributed, therefore we square-root transformed read counts prior to analysis. We ran models for 50,000 iterations, with a burn-in of 5,000 iterations and visually confirmed adequate chain-mixing. Hereafter, we will refer to these estimates simply as consumption.

| Crude protein content of common shrubs
We characterized the nutritional value of common shrubs in each habitat and/or those that were most common in woodrat diets (see below) by measuring relative crude protein content. Crude protein content is considered the best single factor for determining nutritional value of forage plants (Sampson & Jesperson, 1963, pg. 20

| Preference trials
We conducted preference trials in the field from Jun- We calculated a preference index with the following formula:

| Vegetation community
The most common shrubs and trees on the hill were E. nauseosa  Table 1).
Pinus spp. and Phacelia tanacetefolia also occurred in the diet of N.
bryanti with greater than 80% FOO and over 10% RRA in spring and summer combined (Table 1). Neotoma bryanti increased consumption of F. californica in summer relative to spring evidenced by increases in both FOO and RRA (Figure 2, Tables S1 and S2). Neotoma lepida consumed a less diverse diet, with P. fasciculata being the most abundant in spring and summer (FOO = 1.00, RRA = 0.79; Table 1).
Neotoma lepida increased consumption of P. fasciculata from spring to summer (RRA spring = 0.74, RRA summer = 0.91; Figure 2, Tables 1 and S1 and S2). Overall, RRA for the Asteraceae family did not exceed 2% for either N. bryanti or N. lepida and the overall frequency of occurrence was also low (FOO bryanti = 0.37, FOO lepida = 0.12). Thus we are confident that, even with our inability to discriminate within the Asteraceae family, woodrats consume very little if any E. nauseosa or A. tridentata at our site.
Results of our hierarchical Bayesian modeling were consistent with diet composition based on FOO and RRA estimates. Notably, estimates of consumption using bayespref were less extreme than those from average RRA values (Tables 1 and S1-S3, S5, Figure 2).

| Preference trials
A total of 27 individuals were included in diet trials: N. bryanti (n = 15), N. lepida (n = 12). We found that preference was significantly different between species (p < .001,

| D ISCUSS I ON
Despite ongoing hybridization between N. bryanti and N. lepida (Shurtliff et al., 2014), we found differences in dietary preference and dietary composition between these two species; differences that were maintained in both wet and dry years, and across seasons. The primary plants differentially preferred by each species are nutritious relative to other available plants, but also potentially toxic in unique ways, suggesting these species may have evolved or developed distinct metabolic strategies to reduce toxin exposure.
Given the degree of dietary plasticity we observed across seasons in natural diets and in preference trials, we find that N. bryanti is more of a dietary generalist than N. lepida. Dietary differences between the species likely contribute to their spatial segregation across the ecotone, which ultimately determines their opportunities for interspecific interactions, including hybridization.
At this ecotone, the vegetation of the hill community is more diverse and largely distinct from that of the flats, and this diversity and differentiation is partly reflected in the diets of the woodrats that occupy these habitats (Figure 2 and Figure S1, Tables 1 and S1-S3  (Hodgson, 1979).
During spring and summer of 2016, we show that the diet of N.
bryanti on the hill is dominated by F. californica while the diet of N.
lepida on the flats is dominated by P. fasciculata. However, we did find that some N. bryanti on the hill consumed a small amount of P.
fasciculata, while N. lepida on the flats infrequently consumed F. californica. This result from our sample of wild diets is at least partly due to the relative rarity of these two plants in the "alternate" habitat.
However, results of our 2-choice trial show that even when given a choice of both plants, on average, N. bryanti primarily consumed F. californica and N. lepida primarily consumed P. fasciculata. As such, on average, individuals in our experimental trial showed a preference for the plant they most commonly consume in the natural environment. Overall, our field and experimental results demonstrate that N. bryanti show a preference for F. californica and N. lepida show a preference for P. fasciculata, which may reflect differences in behavioral acclimation to different resources and/or underlying species differences in their ability to metabolize these particular plants.
Despite the overall preference N. bryanti and N. lepida exhibit for these plants, there was a great deal of individual variation in our experimental trials. Specifically, most individuals consumed at least some of the presumably novel plant. This is a foraging behavior animals may employ to identify new food resources (Partridge, 1981), and one we might expect when individuals are exposed to novel food items. This short-term consumption of a potentially novel, chemically distinct plant did not appear to have negative consequences for experimental animals as none lost excessive weight over this short period (i.e. >10% body mass) and animals remained alert and respon-   (Jung et al., 2011;Qin et al., 2016). In contrast, P. fasciculata contains cyanogenic glycosides, also highly toxic once cyanide is released from the parent compound (Vetter, 2000). Chemical analysis of plants from the Kelso Valley have shown that F. californica and P. fasciculata contain chemical peaks consistent with the anthraquinone, emodin, in F. californica and the cyanogenic glycoside, prunasin, in P. fasciculata (Matocq et al., 2020). Given the potential toxicity of these plants, why do woodrats eat so much of them? On one hand, from a chemical perspective, these plant species may be among the best of a bad lot. The other common shrubs present-E. nauseosa, A. tridentata, and Ephedra are also known to be chemically well-defended and/or energetically costly to consume (Dial, 1988;Halls et al., 1994;Johnson et al., 1976). In addition to this, though,  (Ulappa et al., 2014).

TA B L E 3 Effects of variables included in linear model of preference trials
The degree to which woodrats and other herbivores can minimize their exposure to toxins by diversifying their diets (Freeland & Janzen, 1974) (Dial, 1988;Shipley et al., 2009;Skopec et al., 2015), and our data support this classifica-  Shurtliff et al., (2013) showed in laboratory trials that the relatively large-bodied N. bryanti is more aggressive than the relatively small-bodied N. lepida. Interspecific competition is thought to be an important driver of dietary differentiation and fine-scale space use in interspecific contact zones between woodrats, with the large-bodied species typically monopolizing optimal nest sites (Cameron, 1971;Dial, 1988). Neotoma bryanti at this site monopolize what is likely the more optimal, relatively thermally stable boulder nesting area of the hill (Brown, 1968). We suspect inherent differences in behavioral, physical, and metabolic capabilities have allowed N. bryanti to monopolize the hill habitat with its diversity of dietary plants, while N. lepida have persisted at the site in part because of its ability to locally specialize on P. fasciculata.
Woodrats are well-known for their capacity to consume large quantities of potentially toxic plants (Larrea tridentata Juniperus sp.-Dial, 1988, Skopec et al., 2007. In particular, N. lepida is known to locally specialize on chemically distinct plants across its range (Larrea tridentata- Mangione et al., 2000; Juniperus- Stones and Hayward (1968), and here, P. fasciculata). The mechanisms that underlie a woodrat's capacity to detoxify these diets likely include expression of their own detoxifying enzymes (Kitanovic et al., 2018;Malenke et al., 2012) and the activity of their gut microbiota (Kohl et al., 2014). Studies are needed to identify loci that are responsible for detoxification of different compounds, the degree to which specific alleles or pathways effectively metabolize particular PSC's, and the interaction between mammalian and microbial genomes in creating toxin resistant phenotypes (Forbey et al., 2018). If unique metabolic adaptations or microbial combinations allow N. bryanti and N. lepida to metabolize different plants in their respective habitats, then migrant and hybrid individuals that do not possess habitat-specific genomic or microbial combinations may suffer reduced fitness (Nosil et al., 2005;Via, 1999;Via et al., 2000).
Selection against migrants would minimize opportunities for interspecific contact and mating (prezygotic isolation), while selection against hybrids with suboptimal allelic or microbial combinations would further limit introgression between the species (postzygotic isolation). Continued integration of field and laboratory studies will be needed to identify the mechanisms that underlie metabolic processing of these diets, and how diet-related selection is influencing the evolutionary trajectory of these species.

ACK N OWLED G M ENTS
We thank Michael Jee, Jason Wurtz, Dan Moore, and Angela Mitchell for assistance in field work. We thank Mozart Fonseca for use of his lab and equipment and Teshome Shenkoru for assistance in measuring plant crude protein content. We thank the members of the UNR Evoldoers lab group for comments on earlier drafts of this manuscript. We thank Mary V. Price and one anonymous reviewer for helpful comments on this manuscript. We also thank M. Forister, J. Forbey, and D. Dearing for helpful discussions. This research was supported in part by the National Science Foundation (IOS-1457209 and OIA-1826801).

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

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw sequences from this study are available in the NCBI Sequence Read Archive under the BioProject accession number PRJNA700875.