Isotopic niche breadth of a generalist mesopredator increases with habitat heterogeneity across its range

. Although generalists are becoming increasingly abundant and widespread, little is known about their response to ecological variation they encounter across their range. For example, the generalist ’ s ﬂ exible diet is cited to help explain recent range expansions, but no study has directly examined this claim. Here, we use stable isotope values of the Virginia opossum ( Didelphis virginiana ), a true generalist, to examine an extension of MacArthur ’ s habitat heterogeneity hypothesis for a single generalist species. If a generalist ’ s diet re ﬂ ects local food abundance, then more heterogeneous landscapes should result in broader niches. We used stable isotope analysis, landcover indices, and WorldClim data to further evaluate how the opossum ’ s use of its environment varies across ancestral regions, expansion fronts, and regions of human-facilitated introductions. Niche breadth varied across its range, especially between expansion fronts. We found a positive relationship between landcover diversity and isotopic niche breadth. WorldClim variables linked to aridity and C 4 plant abundance were most strongly associated with nitrogen ( d 15 N) and carbon ( d 13 C) values, respectively. Our results reveal that a generalist ’ s stable isotope signature re ﬂ ects its local environment, demonstrating their ﬂ exible diet is captured with stable isotopes and supporting the generalist habitat heterogeneity hypothesis.


INTRODUCTION
Scientists have long sought to understand how species are distributed and the underlying forces controlling species ranges. In 1917, the boundaries of an individual species' range were credited to physiological limitations to the environment (Grinnell 1917), and a few years later, Elton defined a species' niche as its biotic position with predator and prey (Elton 1927). To better quantify a species' niche, Hutchinson proposed an n-dimensional hypervolume where the dimensions are the environmental conditions and resources that define a niche (Hutchinson 1957). This method of drawing niches informed the development of analyzing the isotopic niche (Bearhop et al. 2004), typically depicted as an area drawn from two stable isotope ratios that inform a population's environment and diet. Isotopic niches have been leveraged to examine ecological concepts including the niche variation hypothesis and resource breadth hypothesis (Maldonado et al. 2017, Rader et al. 2017. With climate change, human land-use modification, and the extirpation of apex predators, generalists including an array of mammals and birds are becoming more predominant in communities, in some cases expanding their range and entering new communities (Prugh et al. 2009, Davey et al. 2011). Because of their broad diet, generalists can occupy numerous trophic levels within a food web, often making it difficult to predict their impacts on a community or how climate change will impact their populations (Prugh et al. 2009). The literature examining how diet varies across a generalist's geographic range is limited. One study evaluating a generalist species across its range found that variation in diet was often linked to local food abundance (Terraube and Arroyo 2011), while a second study found dietary variation was constricted in highlatitude winters (Soe et al. 2017).
Despite the growing abundance of generalists in communities (Prugh et al. 2009), nothing is known about isotopic niche patterns across a generalist's range. Ongoing range expansions of widespread mesopredators offer natural experiments to investigate isotopic niche variation across a variety of habitats and range locations. The broad geographic range (Fig. 1;Gardner and Sunquist 2003) and northward expansion of the Virginia opossum, Didelphis virginiana (Walsh and Tucker 2018), provide an opportunity to evaluate how a generalist's isotopic niche changes with its environment and at the margins of its range.
The Virginia opossum (henceforth opossum) is a marsupial found from Central America to the Midwestern United States and southeastern Canada. The northward expansion of opossums is well documented in museum records (Walsh and Tucker 2018), and their current northern margins are found west of the Great Lakes in North Dakota (Walsh et al. 2017) and east of the Great Lakes in Maine (Mosby, personal communication). Individual opossums are reported to consume a wide array of foods, that is, the opossum is a type A generalist (Gardner andSunquist 2003, Bearhop et al. 2004).
The habitat heterogeneity hypothesis states that heterogeneous and complex landscapes increase the number of niches and therefore are expected to accommodate more species (MacArthur 1972). Our primary objective was to evaluate a potential extension of the habitat heterogeneity hypothesis for a single generalist species-more heterogeneous landscapes result in the broadening of a generalist species' niche. Our secondary objective was to understand how local climate influences a generalist's isotopic signature.

Sample collection
To examine how the opossum's isotopic niche varies across its range, six regions of interest were defined based on the timing of the opossum's expansion into temperate North America. Mesoamerica and the Gulf Coast were established as ancestral regions based on Pleistocene fossil records (Graham and Lundelius 2010). The Midwest and Northeast were established as recent northward expansion regions, and the California Central Valley and California Coast were established as two regions occupied due to human-facilitated introductions (Fig. 1). Because opossums are relatively nomadic, we conducted our analyses over broad regional spatial scales comparable with the scale required to detect genetic structure (Walsh and Tucker 2018).
Between 20 and 42 individuals were sampled from adult museum skin specimens for each region for a total of 153 specimens sampled (Appendix S1), exceeding the minimum sample size of 5 required per region for stable isotope analyses (Clementz and Koch 2001).

Stable isotope analysis
For each sample, approximately three whole guard hairs were treated with a 2:1 ratio of chloroform-methanol to remove lipids, dried under a fume hood at ambient temperature, diced, and weighed with a Mettler AE 240 balance (Toledo, Ohio, USA) to place 0.5-1.0 mg in a tin capsule. Duplicate capsules for each sample were sent to the University of New Mexico Center for Stable Isotopes (UNM-CSI) for stable isotope analysis of d 13 C and d 15 N. To analyze organic substrates, UNM-CSI uses a Delta V mass spectrometer with a Conflo IV interface (Thermo Scientific, Waltham, Massachusetts, USA), 4010 elemental analyzer (Costech, Valencia, California, USA), and a high-temperature conversion elemental analyzer. Analytical error across 56 runs of the UNM-CSI protein standard (casein) was 0.1099 standard deviations (SD) for d 15 N and 0.0543 SD for d 13 C.
The museum specimens evaluated for this research were collected across 100 yr. Because atmospheric carbon isotope ratios have decreased over the past century due to the influx of greenhouse gases derived from the burning of fossil v www.esajournals.org fuels, we adjusted for time with a linear correction of À0.005& per year between 1917 and 1961 and À0.022& per year after 1961 (i.e., Chamberlain et al. 2005).

Isotopic niche evaluation
The opossums' isotopic niche breadth was calculated for each of the six regions using the Bayesian ellipse method. By implementing Bayesian inference techniques in the R package SIBER v.2.1.3 (Jackson et al. 2011) when calculating the standard ellipse area (SEA B ), sampling error can be accounted for in the isotopic niche.
Because analyses can be partially confounded by sampling regions of unequal geographic area (Stein et al. 2014), two methods were implemented to evaluate the impact that geographic area might have on our calculated SEA B . First, the area of each of the six regions was estimated in ImageJ (https://imagej.nih.gov/ij/) by tracing the polygon created by collection points in a given region (Fig. 1) and was evaluated for correlation with SEA B . For a less coarse approach, SEA B was calculated in SIBER v.2.1.3 (Jackson et al. 2011) for eight subsamples consisting of single or adjoining counties in which at least four opossums were sampled and assessed for correlation with the total land area of the subsamples available from county census records (n = 43; Appendix S1).
To evaluate whether habitat heterogeneity impacts isotopic niche breadth, landcover data from sampled countries were downloaded from https://data.terrapop.org/terraclip#. These landcover data from 2000 assign each geolocated grid cell of 1km resolution one of 22 landcover types (e.g., mixed leaf tree cover, herbaceous cover, broad-leaved evergreen tree cover, mosaic cropland, and tree cover; Fritz et al. 2003). The R script to plot landcover can be found in Data S1. Landcover type for each opossum sample was extracted using its field collection coordinates. The 153 samples were collected from a total of 11 different landcover types (Appendix S1; Appendix S2: Fig. S1). The heterogeneity of the six regions was evaluated using four habitat diversity measurements: Simpson's diversity and evenness indexes and Shannon's diversity and evenness indexes. In our application of Simpson's and Shannon's equations, landcover types were treated as analogous to different species (Lausch and Herzog 2002; Appendix S2: Fig. S1). These diversity measurements (Table 1) were used to evaluate whether habitat heterogeneity of a region is positively correlated with its isotopic niche breadth.

Isotopes and environmental variables
To investigate differences in opossum isotope composition across habitat types, six landcover types were compared (each with n ≥ 7 opossum samples; total n = 133). To investigate the influence that abiotic variables have on opossum isotope composition, 10 climate variables were extracted from the WorldClim dataset (Hijmans et al. 2005) at 10 min resolution using the raster package in R (Hijmans 2015). We used a generalized linear model (GLM) approach to evaluate which climate variables best explained d 13 C and d 15 N patterns. GLMs were designed in IBM SPSS version 26.0 with a normal response for the d 15 N dataset and a gamma response for a positively adjusted d 13 C dataset. Single covariate models were run for each WorldClim variable, and WorldClim variables that were not correlated were also used to build two-covariate and threecovariate models (all variance inflation factors ≤ 1.020). To avoid overfitting and favoring more complex models, model performances were measured with Bayesian information criterion (BIC) and compared by calculating DBIC (BIC i À BIC MINIMUM ; Rafferty 1995, Barker andLink 2015).
In order to compare the opossum, a generalist, to previous stable isotope research on specialist mammals (e.g., Cotton et al. 2016, Smiley et al. 2016, we also evaluated the importance of the ten WorldClim variables for each isotope using conditional forest (CF) analysis with the party package in R (Strobl et al. 2008). CF analysis is a classification and regression tree machine learning method that allows for assessment of non-parametric data and accounts for complex interactions including correlation between independent variables (Strobl et al. 2008). For each isotope, 1000 trees were built using randomly selected WorldClim variables and aggregated. The number of variables to include in each tree was determined based on the forest that yielded the smallest root mean square error measured with the caret R package (Kuhn 2008, Strobl et al. 2008. Because CF analysis does not output

Isotopic niche evaluation
From a sample of 153 opossums collected between 1917 and 2017 (Appendix S1), d 15 N values ranged from 4.35& to 13.98& (AE0.16 SD), and Suess-corrected d 13 C values ranged from À24.99 to À16.23& (AE0.20 SD; Appendix S2: Fig. S2). No region's isotope values were significantly correlated with the year in which they were collected (Spearman's P value ≥ 0.06). The isotopic niche breadth (SEA B ) varied across the opossum's geographic range, with the largest SEA B observed in Mesoamerica and the Midwest, followed by the Gulf Coast (Fig. 2). The smallest SEA B included California's agricultural Central Valley, California's coast, and the Northeast (Fig. 2). There was no correlation between SEA B and county geographic area (Pearson's r = À0.490, P = 0.218; Appendix S2: Table S1). The correlation coefficient increased but remained non-significant when the six regions' SEA B were evaluated for correlation with geographic area (Pearson's r = 0.778, P = 0.069; Table 1). SEA B was positively correlated with Simpson's and Shannon's diversity indexes and Shannon's evenness index (Pearson's r ≥ 0.869, P ≤ 0.025; Table 1).

Isotopes and environmental variables
There were significant differences in d 13 C values between landcover types (Kruskal-Wallis P = 0.021; Appendix S2: Fig. S3). Herbaceous cover had significantly higher d 13 C values compared to broad-leaved evergreen and broadleaved deciduous cover (Dunn's P ≤ 0.029). Similarly, cultivated areas had significantly enriched d 13 C values compared to broad-leaved evergreen and broad-leaved deciduous cover, along with needle-leaved evergreen cover (Dunn's P ≤ 0.046).
For d 15 N, the best performing GLM was mean temperature in the wettest quarter + annual precipitation (BIC = 640.064; parameter estimates = 0.051 AE 0.0186 SE, À0.002 AE 0.0004 SE; P values = 0.006, <0.001 respectively; Table 2). When compared to the best performing model, only one GLM had a DBIC less than two, a value that demonstrates there is only weak support for one model over the other (Rafferty 1995). This competitive model included mean temperature in wettest quarter + precipitation in the driest quarter (DBIC = 0.621; parameter estimates = 0.049 AE 0.0186 SE, À0.006 AE 0.0015 SE; P values = 0.009 and < 0.001 respectively; Table 2). There was strong evidence both of these models For d 13 C, the best performing GLM was maximum temperature (BIC = 593.22; parameter estimate = 0.048 AE 0.0117 SE; P value < 0.001; Table 3). There were one additional GLM with substantial support, which included mean temperature in the wettest quarter (DBIC = 0.684; parameter estimate = 0.021 AE 0.0051 SE; p value < 0.001; Table 3). There was very strong evidence both of these models performed better than the null model (DBIC = 10.917; Table 3; Rafferty 1995).
The conditional forest analysis found that the 10 WorldClim variables selected explained 19.63% of the variance in d 15 N data and 13.53% of the variance in d 13 C data. The WorldClim variable with the highest importance for nitrogen isotopic composition of opossum guard hairs was annual precipitation, followed by precipitation in the driest quarter (Fig. 3A). The WorldClim variable with the highest importance for carbon isotopic composition of opossum guard hairs was precipitation in the coldest quarter, followed by mean temperature in the wettest quarter (Fig. 3B). Precipitation variables were negatively correlated with d 15 N values (Spearman's P ≤ 0.027; Table 4), and temperature variables Notes: P, precipitation; Q, quarter; T, temperature; and BIC, Bayesian information criterion. These models performed significantly better than the null, intercept only, model (Rafferty 1995). Notes: Q, quarter; T, temperature; and BIC, Bayesian information criterion. These models performed significantly better than the null, intercept only, model (Rafferty 1995 (Table 4).

DISCUSSION
We used the opossum, a type A generalist experiencing a range expansion, as a system to study variation in stable isotope values and niche breadth across a broad spatial scale. Temporal expansions and contractions of isotopic niches have been observed from the Holocene through the modern-day samples of small rodents in the Great Basin (Terry 2017). Our results demonstrate that variation in niche size can also occur spatially, as populations span different environmental and climatic conditions. Our analysis supports the hypothesis that the isotopic niche breadth will vary across a generalist's range in correlation with habitat ( Fig. 2; Table 1).
In a study that manipulated barley fields to vary habitat heterogeneity, biologists found that individual diets of generalist arthropods were more variable in the heterogeneous fields. Our range-wide results of a mammalian mesopredator, along with the site-specific study on arthropod generalists (Staudacher et al. 2017), support the generalist habitat heterogeneity hypothesis. To further understand isotopic niche trends in generalists, we suggest sampling type A and type B generalists, as well as specialists, across various degrees of habitat heterogeneity to bolster statistical power and compare results across species. Isotope patterns d 15 N values were primarily driven by precipitation patterns: Areas with lower precipitation and higher temperatures during precipitation extremes had higher d 15 N values (Fig. 3A). Based on our GLM results, for every Celsius degree increase in mean temperature during the wettest quarter, opossum d 15 N increased by 0.051&, and for every millimeter increase of rain, d 15 N decreased by 0.002&. Aridity is known to enrich the nitrogen isotopic composition of soils and vegetation, but this signal can be lost in localized assessment due to other processes that drive d 15 N patterns in microhabitats (Pardo and Nadelhoffer 2010). Our ability to detect a relationship between nitrogen composition and precipitation across a large geographic scale is consistent with a study on two heteromyid rodent species in western North America (Smiley et al. 2016).
Herbaceous and cultivated landcover types, which are often dominated by C 4 plants, were associated with opossums with higher d 13 C values than opossums found in predominantly C 3 forest landcover types (Appendix S2: Fig. S3). Based on our GLM results, for every Celsius degree increase in maximum temperature experienced, opossum d 13 C increased by 0.048&, and for every degree increase in mean temperature during the wettest quarter, d 13 C increased by 0.021&. In general, the WorldClim variables that had the highest importance for d 13 C values (Fig. 3B) and strongest performing models (Table 3) also predict the relative abundance of C 3 and C 4 plant growth-areas with less rain and warmer temperature, climate conditions that are more conducive for C 4 grasses than for C 3 vegetation (Terri and Stowe 1976), had higher d 13 C values. Interestingly, even though opossums are omnivorous generalists, the climate patterns that best predict d 13 C values across their range are consistent with similarly broad spatial studies that evaluated large obligate grazers and small herbivores , Smiley et al. 2016.
WorldClim variables are averaged measurements from 1970 to 2000, and these data are approximate climate variables, especially for v www.esajournals.org opossums collected outside of the 30-yr World-Clim window. The explanatory power of climate variation was low (<20%) for both isotopes in CF analysis, and parameter estimates from the topperforming GLMs were significant but subtle. This may in part be due to sampling from across multiple climatic regimes in North America which could amplify the noise of isotopic signals caused by local ecosystem processes and ecophysiological differences between plant species (Lajtha andMarshall 1994, Pardo andNadelhoffer 2010). This may also be due to the strong influence that human land-use (e.g., urbanization and agriculture) has on mammal communities, especially generalists (Rowe and Terry 2014). The isotopic composition of a consumer is a combination of an animal's foraging behavior and the isoscape in which they are foraging (Yeakel et al. 2016), so it is reasonable for climate not to explain all the variance detected in a relatively nomadic generalist species (Gardner and Sunquist 2003). As a type A generalist (Gardner andSunquist 2003, Bearhop et al. 2004), the opossum's isotopic signature may offer a method to assess the isoscape of broad geographic areas. Our results demonstrate that their d 13 C and d 15 N values are impacted by both the local climate and landcover type. However, because the isotopic values of basal resources vary across space and time (Lajtha and Marshall 1994, Pardo and Nadelhoffer 2010, Yeakel et al. 2016, our dataset cannot be used to directly compare the opossums' diets. To resolve this, primary consumers collected from the same site and year should be analyzed to serve as isotopic baselines to correct each opossum's d 13 C and d 15 N values (e.g., Olsson et al. 2009).

Range expansion
The two regions with expansion fronts, the Midwest and Northeast, did not have similar isotopic niche breadths-the Midwest represented one of the largest niches while the Northeast had the smallest niche (Fig. 2). This difference may be due to the comparative homogeneity of the habitat and climate in the Northeast (Table 1). In addition to encompassing all habitats found in the Northeast, the Midwest included both herbaceous and shrub covers (Appendix S2: Fig. S1). The more heterogeneous habitats in the Midwest (Table 1) that include herbaceous and shrub cover may include a greater variety of C 3 :C 4 plant ratios, resulting in a broader distribution of d 13 C values measured in Midwestern opossums. Similarly, the range of annual precipitation in the Midwest is greater (592 mm range) than in the Northeast (362 mm range). This larger range of precipitation may drive greater variation in d 15 N values in the Midwest. The contrasting isotopic results from the opossum range expansions in the Midwest and Northeast highlight the importance of understanding how local habitat and climate influence the ecology of generalist and mesopredator range expansions.
Because opossums do not have a restricted or preferential diet (Gardner and Sunquist 2003), their stable isotope values should mirror the prey abundance and isotopic variation of the dietary resources in the habitat in which they are foraging. This is congruent with our results-their stable isotopes reflect their local habitat. Further comparative studies could tease apart whether this flexibility is common in both type A and type B generalists and whether this signal is restricted to generalists.

ACKNOWLEDGMENTS
The UNM Center for Stable Isotopes staff analyzed stable isotope ratios. Catherine Badgley, Ben Dantzer, Alison Davis Rabosky, and Tara Smiley provided guidance in design and assessment. This research was supported by a Grant-in-aid of Research from the American Society of Mammalogists, Phi Kappa Phi Society Michigan Chapter, and funds from the University of Michigan's Ecology and Evolutionary Biology Department. We thank the following institutions for providing destructive sample loans of hair (full institution titles in Appendix S1): ASNHC, CCBER, CRCM, CSUC, CSULB, CUMV, FMNH, HSU, LACM, LSU, MCZ, MSB, MSU, MVZ, NYSM, OMNH, OWC, ROM, TTU, UAM, UMMZ, USNM, and UWBM.