Genetic and Environmental Indicators of Climate Change Vulnerability for Desert Bighorn Sheep

Introduction

As evidence for the effects of climate change on biodiversity accumulates (Parmesan and Yohe, 2003; Bellard et al., 2012; Wiens, 2016), assessments of species’ vulnerability to predicted climatic shifts have become increasingly common (Chapman et al., 2014). By providing information on species’ exposure (the magnitude of climate change experienced), sensitivity (the degree to which fitness is affected by a given change in climate), and adaptive capacity (the potential to persist in situ or shift ranges), such assessments can improve our understanding of responses and inform high-level decisions about conserving biodiversity (Dawson et al., 2011). Natural resource agencies are increasingly considering climate vulnerability in management plans (Johnson, 2014; Halofsky et al., 2015; Staudinger et al., 2015), but there are two common characteristics of climate vulnerability assessments that may limit their utility for informing management strategies.

First, climate vulnerability assessments are typically conducted at the species level (Pearson et al., 2014; Pacifici et al., 2015; Urban, 2015), but many of the management decisions that most directly influence species are made at the population level. Natural resource managers must consider how to allocate agency resources among populations within their jurisdiction, and population-level assessments are especially relevant for broadly distributed species of conservation concern because populations are extirpated before species go extinct (Razgour et al., 2018). There is increasing evidence that climate vulnerability can vary considerably among populations (Bay et al., 2018; Razgour et al., 2019) and that failure to consider intraspecific variation can lead to misleading predictions regarding the overall climate vulnerability of a species (D’Amen et al., 2013; Schwalm et al., 2016; Razgour et al., 2019).

Second, vulnerability assessments tend to focus on climate exposure while overlooking other components of vulnerability, such as adaptive capacity (Butt et al., 2016). Assessments often rely on correlative modeling approaches (e.g., ecological niche, bioclimatic envelope, or species distribution models) that use associations between known species occurrences and climatic or other environmental variables to predict the distributions of species across geographic space (Wiens et al., 2009; Anderson, 2013; Porfirio et al., 2014; Foden et al., 2019). These models can be coupled with climate projections to infer climate exposure based on predicted changes in geographic distribution or habitat suitability (Araújo and Peterson, 2012; Foden et al., 2019). However, such assessments typically do not account for the ability of species to adjust to novel climatic conditions via evolutionary adaptation and phenotypic plasticity (Merilä and Hendry, 2014), which can be difficult to measure but potentially critical for understanding how climate change affects species (but see Bush et al., 2016; Razgour et al., 2019). Correlative approaches also rarely consider how effects of climate change on species distribution or persistence may depend on habitat connectivity or changes in landscape configuration (Brooker et al., 2007). These limitations suggest that a more multi-faceted approach to vulnerability assessment may be needed to establish population-level management priorities that better reflect the range of factors influencing vulnerability (Hoban, 2018; Razgour et al., 2018).

Genetic information can inform vulnerability assessments by characterizing the potential for evolutionary adaptation to climate change (Carroll et al., 2014), although whether specific populations can evolve climate-sensitive traits is rarely known (Urban et al., 2016). In the absence of that knowledge, other relevant questions may be addressed with genetic data. (1) How genetically diverse are populations? Genetic diversity is positively correlated with population fitness and persistence for many species (Reed and Frankham, 2003; Frankham, 2005), albeit not universally (Linløkken, 2018), and the amount of standing genetic variation in populations is thought to be the best indicator of their potential for contemporary evolution in response to climate change (Hendry et al., 2011; Sgro et al., 2011). (2) How is genetic variation distributed among populations? Genetically unique populations may be particularly valuable for maintaining species-level evolutionary potential (Petit et al., 2008), and identifying unique populations allows them to be managed as distinct taxa potentially adapted to different local environments (Buchalski et al., 2016). Neutral genetic variation also has been used to distinguish phylogenetic lineages and analyze potential climate-mediated range shifts under the assumption that different lineages likely contain different adaptive diversity (Ikeda et al., 2017), although that approach has been challenged (Smith et al., 2019). (3) How connected are populations via gene flow? Genetic isolation (or conversely, connectivity) of populations may be especially influential for evolutionary adaptation to climate change. Gene flow with neighboring populations can help maintain genetic variation upon which selection acts and can introduce adaptive alleles or gene combinations that confer fitness benefits in novel environments (Sexton et al., 2011). For small, fragmented populations in which effects of genetic drift are strong, gene flow facilitated by connectivity with neighboring populations is likely to be particularly important for both maintaining genetic diversity and allowing spread of alleles conveying fitness advantage (Epps et al., 2005; Creech et al., 2017).

Genetic measures of isolation provide estimates of connectivity and genetic structure, as well as variation in genetic diversity (Hedgecock et al., 2007; Lowe and Allendorf, 2010; Meirmans and Hedrick, 2011). However, genetic measures of isolation may become decoupled from landscape-based measures of isolation as a result of wildlife translocations (Malaney et al., 2015; Jahner et al., 2019) or rapid landscape change (Epps and Keyghobadi, 2015). Consequently, genetic-based isolation measures may poorly reflect the degree to which populations are linked via natural dispersal and subsequent gene flow with neighboring populations in the contemporary landscape. In such cases, spatial measures of population isolation may be more useful, such as landscape resistance modeling (Epps et al., 2007; Chetkiewicz and Boyce, 2009) or even Euclidean distance among habitat patches (Prugh, 2009).

As a case study, we illustrate how genetic and environmental data can be combined to assess population-level vulnerability and inform management priorities for desert bighorn sheep (Ovis canadensis nelsoni) within and surrounding lands administered by the U.S. National Park Service (NPS). Desert bighorn sheep, one of three subspecies of North American bighorn sheep (Wehausen and Ramey, 2000), range across the southwestern U.S. and northern Mexico. Bighorn are habitat specialists that rely on steeply sloped, rugged terrain with sparse vegetation and good visibility to detect and evade predators (Gionfriddo and Krausman, 1986; Sappington et al., 2007). Much of the best remaining habitat and many of the largest populations of desert bighorn are located on federal lands, including lands managed by NPS (Figure 1), and these federal lands have greater exposure to climate change than the U.S. as a whole (Gonzalez et al., 2018). Bighorn population dynamics are closely tied to temperature and precipitation patterns (Wehausen et al., 1987; Douglas, 2001; Epps et al., 2004; Bender and Weisenberger, 2005). Projected changes to the climate of the southwestern U.S. during this century include higher annual average temperatures, longer and hotter summer heat waves, reduced winter and spring precipitation, less reliable surface water supply, and more frequent and severe droughts (Garfin et al., 2014). If manifested, such climatic changes could negatively affect bighorn sheep populations by increasing physiological stress or by reducing availability of forage or water. Yet, evaluating the effect of these changes using a correlative approach would be challenging because of the definitive influence of landscape configuration and isolation on genetic diversity and gene flow of bighorn populations (Epps et al., 2005, 2006; Creech et al., 2017). Moreover, potential climate-related adaptive variation among desert bighorn sheep populations has been identified (Buchalski et al., unpublished) that suggests responses may differ in different habitats.

We developed an extensive genetic dataset for desert bighorn sheep across 10 NPS units and surrounding lands, and we used population-level measures of genetic diversity and region-wide analyses of genetic structure and gene flow to estimate relative capacity of populations to cope with climate change via evolutionary adaptation. We complemented genetic analyses with estimates of population isolation based on empirically derived landscape resistance models, which may better characterize isolation of desert bighorn populations with a history of translocation. We estimated climate exposure of populations using a species-neutral climate change metric (forward climate velocity; Loarie et al., 2009) and an environmental variable associated with persistence probability of desert bighorn populations (maximum elevation of habitat patch; Epps et al., 2004).

We took advantage of the history of population augmentation and reintroduction efforts for desert bighorn sheep in portions of our study area to explore the influence of translocations on genetic diversity and genetic structure. We expected that translocated populations would not be in genetic equilibrium with nearby populations because their genetic makeup would reflect influences of the artificial movements of individuals from potentially distant populations, rather than (or in addition to) influences of natural gene flow with neighboring populations constrained by landscape characteristics and time. Because of this expected decoupling of geographic influences from gene flow and resulting genetic diversity in translocated populations, we predicted that: (1) populations that were more geographically isolated (i.e., further from neighboring populations) would have lower genetic diversity, but this negative relationship would be stronger for native than for translocated populations; (2) populations that were more genetically isolated (i.e., more genetically differentiated from nearby populations) would have lower genetic diversity, but this negative relationship would be stronger for native than for translocated populations; and (3) populations that were more geographically isolated would be more genetically isolated, but this positive relationship would be stronger for native than for translocated populations.

Finally, we considered the implications of our results for management of desert bighorn sheep, including how the information from this study could contribute to assessments of the relative vulnerability of desert bighorn populations to climate change effects.

Materials and Methods

Study Area

The study area encompassed 10 NPS units and adjacent lands (e.g., state lands, Bureau of Land Management lands, U.S. Forest Service lands, or Indian reservations) that contained bighorn populations likely to interact with those on NPS lands (Figure 1). This study area included the majority of desert bighorn sheep on NPS lands and was heterogeneous with respect to many landscape characteristics that influence bighorn sheep. In all areas, bighorn habitat was defined by steep terrain and proximity to reliable surface water. However, the configuration of such habitat varied considerably, including areas where habitat was very discretely distributed (e.g., the Mojave Desert) and areas where habitat was relatively continuous (e.g., the Grand Canyon area). Three arid regions with different climate regimes and biota were represented: the Mojave and Sonoran deserts and the Colorado Plateau (Bender, 1982). In most areas the landscape was minimally altered by anthropogenic development (urbanization, highways, mining, water impoundments, etc.), but the extent of developments was greatest in southern California and Nevada. The degree to which population history of bighorn sheep had been directly influenced by management actions also varied. Populations in California, Arizona, and southern Nevada were predominately native, while Utah contained many populations that were reintroduced during the past half century using individuals sourced from distant areas in some cases (Epps et al., 2003; Utah Division of Wildlife Resources, 2013; Wild Sheep Working Group, 2015; Jahner et al., 2019).

Genetic Sampling and Genotyping Methods

We primarily used non-invasive sampling of fecal pellets to obtain DNA from individuals (Wehausen et al., 2004) and combined genetic datasets from multiple projects covering different portions of the study area. Populations were sampled during 2000–2003 in the southern Mojave Desert (Epps et al., 2006), 2003–2007 in southern Nevada and near Lake Mead (Jaeger and Wehausen, 2012), and 2003–2013 in Death Valley National Park and surrounding habitat. These sampling efforts generally targeted waterholes where bighorn sheep congregated during the summer months. Populations in northern Arizona and southern Utah were sampled during 2011–2014 with survey, radiotelemetry, and sightings data from state wildlife agencies and NPS used to identify areas in which to focus efforts. DNA sources also included a small number of blood and tissue samples (<5% of all samples) collected from carcasses discovered in the field or contributed by managing agencies, including California Department of Fish and Wildlife, Nevada Department of Wildlife, Arizona Game and Fish Department, U.S. Geological Survey, and NPS. Differences in sampling periods among populations were unlikely to significantly influence our findings because (1) dynamics of genetic structure and diversity operate on a generational basis, and the 14-year period between the earliest and latest sampling efforts represented only approximately two bighorn generations, (2) no major changes to the landscape (e.g., highway or dam construction) occurred during this period that would strongly impact gene flow, and (3) follow-up research on a subset of Mojave Desert populations showed no significant changes in genetic diversity over the course of a ∼12-year period (Epps et al., 2018).

Samples were processed and genotyped in three different laboratories (corresponding to the different projects) using similar techniques. Each lab genotyped samples at a slightly different set of neutral or neutral-acting microsatellite loci (14–17 loci per sample, depending on the laboratory, with 10 loci common to all laboratories). Nickerson (2014) tested both linked and putatively neutral loci for evidence of selection in bighorn populations within the Mojave Desert region using outlier tests, and only locus BL4 (linked to interferon gamma gene) showed potential evidence of positive selection. Evaluation of later BL4 data in the Mojave region (Epps et al., 2018) showed that the signal of positive selection for BL4 was transitory (i.e., no longer present two bighorn sheep generations later). Due to the high strength of genetic drift in these small populations (Creech et al., 2017), even loci linked to known genes largely reflect patterns at putatively neutral loci (Nickerson, 2014). To realign allele sizes for consistency across laboratories, we genotyped a small subset of samples at each locus used by multiple laboratories and translated all allele sizes to match those used in the most recent project. A full description of genetic laboratory protocols for the most recent project (samples from 2011 to 2014) can be found in Supplementary Material, and protocols for earlier projects can be found in Epps et al. (2006) and Jaeger and Wehausen (2012).

Defining Populations and Habitat Polygons

We grouped individuals into populations based on the location where they were sampled, then created an associated habitat polygon representing the geographic area inhabited by each population for use in spatially explicit analyses. In the Mojave and Sonoran deserts (southern California, southern Nevada, and southwestern Arizona) habitat was generally distributed in discrete patches of steeply sloped terrain separated by broad valleys, which made assigning individuals to populations and mapping associated habitat polygons relatively straightforward. We used a 10% slope cutoff to establish the boundaries of habitat polygons in these regions (Epps et al., 2007) and relied on expert opinion to modify polygon boundaries in areas where this cutoff did not accurately represent the extent of habitat known to be used by a population.

In the Colorado Plateau region (northern Arizona, southern and southeastern Utah) habitat was more continuously distributed and establishing populations and habitat polygon boundaries was less straightforward. We defined populations based on the spatial clustering of individuals (i.e., groups of samples whose locations were clearly separated from other groups) and used the spatially explicit genetic clustering program GENELAND (Guillot et al., 2005) to assign individuals to populations when we were unsure of population boundaries. After establishing populations, we created habitat polygons by generating a minimum convex polygon (MCP) from sample locations, adding a 10 km buffer to reflect the likely use of areas beyond the MCP, and removing areas inside the buffered MCP that were not suitable habitat (e.g., flat terrain, water bodies).

We used bighorn sheep translocation records from the Western Association of Fish and Wildlife Agencies (Wild Sheep Working Group, 2015) to determine whether populations were native or had a history of translocation (reintroduction or augmentation with individuals from other populations) to aid in interpretation of genetic structure results.

Genetic Structure

Individuals were distributed relatively continuously rather than in discrete populations in portions of our study area, and preliminary analyses suggested a pattern of genetic isolation by distance (Supplementary Figure S1) that could produce misleading patterns of genetic structure (Kalinowski, 2011; Janes et al., 2017; Perez et al., 2018) when analyzed using common Bayesian clustering algorithms such as STRUCTURE (Pritchard et al., 2000). We therefore relied on two analytical methods expected to perform more reliably when applied to genetic datasets that deviate from the classic island model of population structure. First, we used discriminant analysis of principal components (DAPC), a multivariate method that summarizes between-group genetic differentiation while ignoring within-group variation and making no assumptions about the underlying population genetic model (Jombart et al., 2010). We conducted DAPC using the adegenet package (Jombart, 2008) for R (R Development Core Team, 2017) and determined the most likely number of genetic clusters using 20 replicate runs of K-means clustering. We used alpha-score optimization to determine the number of retained principal components that represented the best trade-off between discrimination power and overfitting. We used a scatterplot of the first two discriminant functions to assess relatedness between clusters and estimated population-level cluster assignment probabilities by averaging individual assignment probabilities for the most likely number of clusters. Initial results suggested a major genetic split between one cluster and all remaining clusters, and we repeated the analysis after excluding the divergent cluster to identify any additional hierarchical structure.

We used the clustering program TESS3 (Caye et al., 2018) as an alternative approach to determine genetic structure. This program requires no assumptions about linkage or Hardy-Weinberg equilibrium and can incorporate information on sample locations using a spatial regularization parameter λ that adjusts the strength of spatial dependence (i.e., the degree to which geographically proximate individuals are more likely than distant individuals to share ancestral genotypes). We considered λ values spanning a range of spatial dependence strengths, allowing the number of clusters (K) to vary from 1 to 20, with 50 replicates for each K-value. We retained the run with the lowest root mean squared error for each K-value and examined cross-validation plots based on a cross-entropy criterion (Frichot and François, 2015) to determine the optimal number of clusters for each λ value tested. All genetic structure analyses were run using a subset of our genetic data that included 10 loci used by all labs and for all geographic regions, thus avoiding potential bias in these analyses caused by missing data.

To confirm the pattern of isolation by distance among populations in our study area, we calculated geographic distances (Euclidean distances between centroids of populations’ habitat polygons) and genetic distances (Weir and Cockerham (1984) F_ST) for all pairs of populations. We then performed a Mantel test of the relationship between log-transformed geographic distance and linearized genetic distance (F_ST/[1-F_ST]; Rousset, 1997).

Genetic Diversity

We assessed genetic diversity using genotypes at neutral microsatellite loci. For each population we calculated allelic richness (A_r), the average number of alleles per locus after correcting for variation in sample sizes among populations, using rarefaction (Leberg, 2002) with a minimum sample size of six individuals. To facilitate comparisons with populations from other studies, we also calculated expected heterozygosity (H_e), a common genetic diversity metric that is insensitive to sample size. We used the gstudio package (Dyer, 2014) for R to calculate both metrics.

Isolation

We quantified the isolation of each population using both geographic (i.e., landscape) and genetic measures. We used both approaches because genetic differentiation metrics based on direct measures of gene flow among populations such as F_ST are strongly influenced by time lags, translocations, and other non-equilibrium situations. In particular, the extensive history of translocations in parts of the study area made F_ST a questionable indicator of natural gene flow among some populations, and thus of limited value for predicting how landscape structure will influence future population connectivity. We assessed geographic isolation using a landscape resistance-based approach rather than Euclidean distance to account for the fact that populations geographically close to each other but separated by terrain that is highly resistant to dispersal (e.g., broad valleys or highways) may be less connected than geographically distant populations separated by terrain that is favorable for dispersal. We used a resistance model previously developed for desert bighorn sheep in our study area that included a Gaussian effect of slope (where moderate slopes were less resistant than shallow or steep slopes) and strong barrier effects of major water bodies and interstate highways (Creech et al., 2017; Supplementary Material). Terrain ruggedness is also known to influence bighorn habitat use (Sappington et al., 2007) and potentially dispersal but was not included in the resistance model because it was very highly correlated with slope at the spatial resolution (∼100 m) at which the resistance model was developed. We used the model to calculate effective distance (a measure that combines geographic distance and relative habitat resistance) along the least-cost path between the borders of each pair of habitat polygons using Linkage Mapper version 1.1.1 (McRae and Kavanagh, 2011). We quantified geographic isolation of each population as the harmonic mean of effective distances to the nearest three neighboring populations (Epps et al., 2005). We then quantified genetic isolation of each population as the mean pairwise F_ST (Weir and Cockerham, 1984) with these three nearest neighboring populations determined by effective distances.

Relationships Among Genetic Diversity and Isolation Measures

We used simple linear regression to test our predictions that (1) native populations would show stronger negative relationships between genetic diversity and geographic isolation than translocated populations, and (2) native populations would show stronger negative relationships between genetic diversity and genetic isolation than translocated populations. For each combination of genetic diversity (allelic richness) and isolation measure (mean pairwise F_ST or mean effective distance to three nearest neighboring populations), we fit separate regression models for native (n = 48) and translocated (n = 14) populations. We also fit separate regression models of the relationship between genetic isolation and geographic isolation for native and translocated populations to determine whether the degree of concordance between these isolation measures was affected by translocations. We compared model r² and P-values to determine the relative strength of relationships for the two types of populations.

Climate Change Exposure

We used forward climate velocity as an index of the relative exposure of bighorn populations to climate change (Carroll et al., 2015). Forward climate velocity is a measure of the local direction and speed of climate change that indicates how species would need to shift distributions to remain within climatic envelopes (Dobrowski et al., 2013; Hamann et al., 2015). It is a species-neutral measure of exposure that may serve as a useful alternative to more complex correlative models when climate variables influencing the geographic distribution of a species are poorly understood or are believed to vary across the range (Brito-Morales et al., 2018). We obtained gridded spatial data on forward climate velocity based on mean projections for the 2050s from an ensemble of 15 CMIP5 global climate models (AdaptWest Project, 2015). We considered two emissions scenarios: (1) RCP4.5, a moderate emissions scenario that assumes climate policies are able to stabilize climate change by 2100, and (2) RCP8.5, a high emissions scenario that assumes that greenhouse gas emissions continue to rise throughout the twenty-first century (van Vuuren et al., 2011). Climate velocities were derived from a principal component analysis of 11 biologically relevant climate variables, and thus represent predicted potential shifts in overall climate rather than shifts in any single climate variable. We believe this approach was suitable for assessing relative exposure of desert bighorn populations across a large geographic range because although temperature and precipitation are known to influence fitness of desert bighorn sheep, the specific climatic conditions to which bighorn are most sensitive are not fully understood and may vary geographically. For each emissions scenario, we calculated the mean forward climate velocity for grid cells within each population habitat polygon as an estimate of predicted exposure to climate change.

Additionally, we calculated the maximum elevation within each habitat polygon from a 10 m resolution digital elevation model as an indicator of the relative climate stress currently experienced by desert bighorn populations. Previous research in the Mojave Desert region of California found higher extinction risk and lower genetic diversity in populations occupying mountain ranges with lower maximum elevation, which is typically associated with higher temperature and less precipitation (Epps et al., 2004, 2006). We used simple linear regression to explore whether any relationship existed between patch-level estimates of the two measures of climate change exposure (forward climate velocity and maximum elevation).

Overall Vulnerability

To aid in interpretation of results across multiple variables related to climate vulnerability, we created an index of overall vulnerability. We converted population values for genetic diversity (using A_r values only), genetic isolation, geographic isolation, forward climate velocity (using RCP4.5 values only), and maximum elevation to percentiles, reversing the order of values as necessary to ensure that higher percentiles represented greater vulnerability for all variables. We then calculated the mean percentile across the five variables for each population as an index of that population’s overall vulnerability with respect to the components of vulnerability included in our analysis. We recognize that these components likely vary in their relative influence on overall vulnerability, but we assigned them equal weight as a potentially useful starting point for making comparisons among populations.

Results

We genotyped 1,652 individuals from 62 populations (Figure 1). DAPC results indicated that the most likely number of clusters varied from four to 28, with a modal value of 11 across replicate runs (Supplementary Figure S2), suggesting that the program had difficulty reliably determining the number of genetic clusters. We generated individual- and population-level assignment probabilities assuming K = 11, and results revealed a single cluster corresponding to populations in southeastern Utah that was differentiated from all other clusters further west (Figure 2, upper panel). Nine populations had highest assignment to this distinct cluster, including seven of the 14 translocated populations within our study area. When these nine populations comprising the distinct cluster were excluded from the dataset, reanalysis suggested that remaining populations comprised 10 clusters that were not highly distinct, with most populations showing strongly mixed assignment (Figure 2, lower panel).

TESS3 results suggested that the most likely number of clusters was five, regardless of the strength of spatial dependence assumed by the model (Figure 3 and Supplementary Figure S3). Clusters roughly corresponded to the following geographic groupings: (1) southeastern Utah populations within and around Arches National Park (ARCH), Canyonlands National Park (CANY), Capitol Reef National Park (CARE), and Glen Canyon National Recreation Area (GLCA); (2) southern Nevada/western Arizona populations within and around Lake Mead National Recreation Area (LAKE), along with Zion National Park (ZION); (3) northern Arizona populations within and around Grand Canyon National Park (GRCA); (4) southern California populations north of Interstate 15, within and around Death Valley National Park (DEVA); and (5) southern California populations south of I-15, within and around Mojave National Preserve (MOJA) and Joshua Tree National Park (JOTR). Ancestry was mixed for most populations, with assignment probabilities of > 10 percent to two or more clusters (Figure 3), and the vast majority of individuals within populations also exhibited assignment probabilities of > 10% to multiple clusters (Supplementary Figure S4). Most populations with highest assignment to the southeastern Utah cluster, and nearly half of populations with highest assignment to the southern Nevada/western Arizona cluster, were translocated. Populations with highest assignment to the remaining clusters were almost exclusively native. A Mantel test confirmed that populations exhibited genetic isolation by distance across the study area (r = 0.21, p < 0.001).

Allelic richness across the 62 populations ranged from 2.32 to 3.90 with a mean of 3.24, while expected heterozygosity ranged from 0.44 to 0.70 with a mean of 0.60 (Table 1). A_r and H_e of populations were highly correlated (r = 0.95). We focus henceforth on A_r results because A_r is sensitive to population bottlenecks and is considered a better indicator of evolutionary potential than other genetic diversity metrics (Leberg, 2002; Allendorf and Luikart, 2009). We observed the highest genetic diversity in populations located within or near eastern DEVA, LAKE, and GRCA (Figure 4A). Populations in southeastern Utah exhibited the lowest genetic diversity (especially the Arches-Gemini Bridges and Island in the Sky populations), as did many populations in the southern Mojave Desert (e.g., N. San Bernardino Mts./Cushenbury Mts., Newberry Mts./Ord Mts./Rodman Mts., and Iron Mts. populations).

Mean pairwise F_ST with neighboring populations varied from 0.03 to 0.30 (Table 1). The most genetically isolated populations (i.e., those with highest pairwise F_ST with neighboring populations) occupied island ranges in the Mojave Desert, including the Newberry Mts./Ord Mts./Rodman Mts., N. San Bernardino Mts./Cushenbury Mts., San Gabriel Mts., and Turtle Mts. populations (Figure 4B). Populations near the California-Nevada border in the area between MOJA, DEVA, and LAKE (e.g., Eldorado Mts., N. Spring Range, Clark Mts./S. Spring Range populations) were the least genetically isolated from their neighbors.

The most geographically isolated populations (i.e., those with highest effective distances to neighboring populations) were located in the southern Mojave Desert, primarily outside of MOJA and JOTR (e.g., Chemehuevi Mts., Orocopia Mts., Cady Mts. populations), as well as the Zion population in southern Utah (Figure 4C). These populations tended to occupy island ranges surrounded by desert flats that were unfavorable to dispersal. In contrast, populations within and around DEVA, GRCA, and southeastern Utah, where steeply sloped terrain was more continuous, were relatively well connected geographically to their neighbors.

Genetic diversity of populations was negatively associated with degree of genetic isolation (Figure 5), and this relationship was stronger for native (r² = 0.53, p < 0.001) than for translocated (r² = 0.29, p = 0.047) populations. Genetic diversity was also negatively associated with geographic isolation (Figure 6) for native populations (r² = 0.27, p < 0.001), however, there was no relationship between genetic diversity and geographic isolation for translocated populations (r² = 0.007, p = 0.769). Geographic and genetic measures of isolation were positively associated (Figure 7) for native populations (r² = 0.44, p < 0.001), but this association was weak for translocated populations (r² = 0.12, p = 0.220).

Mean forward climate velocities of populations were on average 50% higher for the RCP8.5 emissions scenario than for the RCP4.5 scenario, but velocities for the two scenarios were highly correlated (r = 0.95; Supplementary Figure S5). Because we were primarily interested in the relative differences in velocity among populations, we henceforth focus on results for the RCP4.5 scenario (but see Supplementary Figure S6). The highest forward climate velocities were observed for the southernmost populations in the study area within and around JOTR, particularly the Eagle Mts. and Orocopia Mts. populations, and for the White Mts. and Inyo Mts. populations in the northern Mojave Desert (Figure 4D). The lowest forward climate velocities were observed for populations within and around GRCA and LAKE, with low-to-moderate values for populations in the northeastern portion of GLCA, ZION, portions of DEVA, and the San Gabriel Mts. population in the western Mojave Desert.

Maximum elevation of habitat polygons ranged from 987 to 3,498 m (Table 1). Populations in the southern Mojave Desert from JOTR and MOJA east to the Colorado River (e.g., Iron Mts., S. Bristol Mts., Chemehuevi Mts. populations) had the lowest maximum elevations, while populations in the northern Mojave Desert (e.g., White Mts. and Inyo Mts. populations) and western Mojave Desert (e.g., San Gabriel Mts. and San Gorgonio Mts. populations) had the highest maximum elevations (Figure 4E). Maximum elevation and mean forward climate velocity were not related to each other (r² = 0.003, p = 0.64), although a weak negative relationship was observed if one extreme outlier population (White Mts.) was removed from the analysis (r² = 0.12, p = 0.006; Figure 8).

When vulnerability results were combined across all variables (i.e., genetic diversity, genetic isolation, geographic isolation, climate velocity, and maximum elevation), the populations with highest overall vulnerability were primarily located in the southern Mojave Desert and in southeastern Utah. Populations in and around DEVA and GRCA exhibited the lowest overall vulnerability (Figures 4F, 9).

Discussion

Much emphasis has recently been placed on identifying and prioritizing species that are most vulnerable to climate change (Thomas et al., 2004; Williams et al., 2008; Summers et al., 2012; Foden et al., 2013). In contrast, our analysis yielded information that may be useful in prioritizing conservation actions for populations of a single species, an approach that has been employed much less frequently in the context of climate change (but see Blair et al., 2012; Razgour et al., 2018). We found wide variation among desert bighorn sheep populations with respect to several key genetic and environmental characteristics expected to influence their vulnerability to climate change. Our analysis showed local and regional correspondence along several axes of vulnerability: for instance, the southernmost populations in the study area had the highest forward climate velocities; in that region, likewise, elevations were lowest and many populations had high isolation and low genetic diversity (Figure 4). Other areas showed low vulnerability due to low forward climate velocity, high elevation, and high genetic and geographic connectivity, such as DEVA and GRCA (Figure 4). These assessments provide important context for management of desert bighorn sheep across the region. By incorporating estimates of genetic structure and diversity, gene flow, and landscape connectivity, we were able to characterize evolutionary aspects of climate vulnerability that are often missing from commonly used correlative approaches.

Our assessment of genetic structure incorporated nearly all native populations of desert bighorn sheep found within the “Nelson” lineage described by Buchalski et al. (2016), which encompasses populations in the Mojave Desert and parts of the Great Basin and Colorado Plateau ecoregions. Concordant with being part of a single lineage of desert bighorn sheep, we found that populations throughout our study area exhibited mixed ancestry, with nearly all populations (and individuals within populations) assigning to multiple genetic clusters with non-negligible probabilities in both TESS3 and DAPC analyses (Figures 2, 3). This is consistent with the pattern of isolation by distance observed in our data set. Previous studies have in some cases suggested clearer genetic divisions, but those studies were conducted at smaller spatial extents or used different analytical approaches aimed at distinguishing recent gene flow. For instance, Jaeger and Wehausen (2012) found almost no gene flow between populations on opposite sides of the Colorado River in the LAKE region; yet our analysis suggested that populations on opposite sides of the Colorado River within either the LAKE or GRCA region were more genetically similar than were geographically distant populations located on the same side of the river. In another case, the strong genetic separation described by Epps et al. (2018) between populations north and south of Interstate 40 in southern California was much weaker according to our analysis (Figures 2, 3). These conflicting results are suggestive of hierarchical genetic structure. The previous studies conducted at smaller spatial extents revealed lower-level structure, while our analysis, which included more of the total genetic variation present within the subspecies, revealed higher-level structure. The seemingly contradictory findings may also speak to the relatively recent influence of dams on gene flow across and along the Colorado River. Prior to construction of dams in the early-to-mid twentieth century, bighorn sheep would have been able to cross the Colorado River relatively easily during periods of low streamflow. The large reservoirs and dam-regulated streamflow in recent decades appear to have largely prevented such crossings and increased genetic divergence of populations on opposite sides of the river. Thus, our current large-scale analysis reflecting longer-term patterns did not indicate a strong barrier effect of the Colorado River, while analyses using methods that characterize recent gene flow (e.g., assignment tests) suggest otherwise (Creech et al., unpublished; Wehausen et al., unpublished).

The most significant genetic division in our study area appeared to be between populations in southeastern Utah (in and around ARCH, CANY, and CARE) and all other populations (Figure 2). The history of the populations in this region provides a plausible explanation for this pattern. Nearly all populations in southeastern Utah went extinct between the late 1800s and the 1940s as a result of livestock-borne diseases, unregulated harvest, and mining activities (Singer and Gudorf, 1999), leaving only a small remnant population in CANY. This population bottleneck likely reduced genetic diversity and interrupted gene flow with adjacent regions, causing allele frequencies to diverge, and resulting in the genetic distinction of this remnant population. Reintroductions of bighorn sheep populations to southeastern Utah beginning in the mid-1970s using individuals sourced from the remnant CANY population (Utah Division of Wildlife Resources, 2013) likely did little to change the genetic distinction of this region. Other signatures of past translocation events were also present in the genetic structure results; for instance, the Zion population was most similar to populations in LAKE because it was reestablished using individuals translocated from the River Mountains in Nevada. Although southeastern Utah populations exhibited relatively low isolation and low forward climate velocity, the loss of genetic diversity due to demographic history and reintroduction significantly elevated their overall vulnerability in our assessment (Figure 4).

Genetic diversity varied strongly among populations, but was high throughout much of our study area in comparison to estimates from other parts of the range of bighorn sheep, suggesting that potential for evolutionary adaptation in these mostly native populations remains high despite small population sizes and rapid genetic drift (e.g., Epps et al., 2005). GRCA, DEVA, and LAKE contained some of the most genetically diverse bighorn sheep populations reported in the literature. Populations occupying these NPS units had mean expected heterozygosities of 0.64, 0.63, and 0.65, respectively. These estimates are higher than microsatellite-based H_e estimates for most desert bighorn sheep populations outside our study area: Boyce et al. (1997) reported mean H_e of 0.55 for populations in the Peninsular Ranges of southern California and 0.50 for populations in southern New Mexico; and Gutiérrez-Espeleta et al. (2000) found mean H_e of 0.57 for populations in southern Arizona. H_e of populations in our study area also compared favorably to reported mean H_e for Rocky Mountain bighorn sheep (O. c. canadensis) populations of 0.57 (Gutiérrez-Espeleta et al., 2000), 0.55 (Forbes et al., 1995), and 0.63 (Driscoll et al., 2015). We note, however, that these previous studies used sets of microsatellite loci that were similar but not identical to those used in our study, which could potentially influence these comparisons.

Our results also suggest that genetic diversity of some populations in our study area may be in flux because translocated populations have not yet reached equilibrium. As predicted, geographic isolation was negatively correlated with genetic diversity for native but not for translocated populations (Figure 6), and the negative correlation between genetic isolation and genetic diversity that we observed was stronger for native than for translocated populations (Figure 5). The stronger correlation between genetic isolation and geographic isolation for native than for translocated populations (Figure 7) was also consistent with our predictions. These results suggest that despite the small census sizes and rapid genetic drift typical for desert bighorn sheep populations (Epps et al., 2003, 2005), translocated populations do not appear to have achieved equilibrium with respect to influences of landscape connectivity on genetic diversity and genetic differentiation. Where translocation events have muddied the relationship between landscape characteristics and gene flow, geographic measures of isolation may provide a more accurate picture of potential for future gene flow. Geographic isolation measures may also be more useful where populations have not yet reached genetic equilibrium following major changes to the landscape (e.g., highway or dam construction). We therefore suggest that managers consider both geographic and genetic measures when assessing population isolation in regions where translocation efforts have occurred.

We assessed genetic diversity of populations using neutral microsatellite markers because a major goal of our research was to explore landscape influences on gene flow, and neutral markers provide unbiased estimates of demographic processes, such as gene flow and genetic drift (Holderegger et al., 2006). Evolutionary potential, however, depends on adaptive genetic variation – variation at genes that affect fitness – rather than neutral variation. While the correlation between neutral genetic variation and quantitative variation cannot always be assumed (Reed and Frankham, 2001), genetic diversity at neutral and adaptive-linked microsatellite loci were shown to be strongly correlated for desert bighorn sheep populations in the Mojave Desert (Nickerson, 2014). Heterozygosity estimates for bighorn sheep derived from microsatellites and SNPs were also correlated (Miller et al., 2014). Thus, we believe that our genetic diversity estimates may be useful if imperfect indicators of evolutionary potential of desert bighorn sheep populations. Genomic methods such as outlier tests and genotype-environment association analyses (e.g., Prunier et al., 2011; Razgour et al., 2018) could provide more direct estimates of adaptive genetic variation and also shed light on relative sensitivity of populations to climate change by revealing populations harboring alleles associated with higher fitness in novel climates. Yet, one such recent effort determined that the strong hierarchical spatial structure observed in bighorn sheep confounds interpretation of adaptive differences across many habitats (Buchalski et al., unpublished).

The forward climate velocity estimates in our analysis are subject to considerable uncertainty associated with climate model projections and emissions scenarios. This issue affects all vulnerability analyses that rely on climate projections. However, several aspects of our analysis give us confidence in our velocity-based conclusions regarding climate exposure. First, the velocity data we used were generated using a diverse set of climate variables and an ensemble of climate models, minimizing the effects of individual variables or models that behave as outliers. Second, we considered mid-term (2050s) climate projections, avoiding the increased uncertainty associated with longer-term projections. Third, we considered two emissions scenarios that bracket a wide range of possible climate trajectories. Fourth, our conclusions were based on relative differences in velocity among populations, which were largely consistent between emissions scenarios, rather than absolute magnitudes of climate velocity.

The weak relationship observed between our two measures of climate change exposure is not necessarily surprising because they were chosen to represent different components of exposure: maximum elevation as an index of current climate stress, and forward climate velocity as an index of the magnitude of expected climatic change. In some cases, these metrics indicated opposite conclusions regarding relative exposure of populations, as illustrated by the White Mts. population near the northwestern limit of our study area. This population had the highest maximum elevation, indicating low climate stress at present because the population presumably experiences relatively cool and wet conditions. The White Mts. population also had the highest mean forward climate velocity, indicating high exposure over the long term, because it already occupies the highest-elevation habitat in the vicinity and long-distance dispersal would be necessary to find suitable future habitat at even higher elevation or higher latitude. Interestingly, a recent analysis of adaptive genetic variation in desert bighorn sheep suggests that the White Mts. population has a high frequency of an allele associated with high-elevation habitat (Buchalski et al., unpublished); high-elevation habitat presents different life-history challenges due to altitude and cold stress, suggesting that local adaptation could be expected. The question of whether effects on populations will be greater in areas experiencing the largest climatic changes or in areas where climatic conditions are closest to tolerance limits has not been well explored (but see Beever et al., 2010), despite having important implications for management.

We urge caution in interpreting maximum elevation results at a range-wide scale. Evidence for the relationship between maximum elevation and population persistence comes from populations in the southern Mojave Desert, where elevation is strongly associated with precipitation and forage availability (Epps et al., 2004). Whether elevation is as strong a determinant of population persistence in other regions is not clear, nor is it clear whether cross-regional comparisons are appropriate given major differences in climate, vegetation type, and topography. We suspect that maximum elevation is most useful as a measure of differences in climate change exposure among populations within a region.

Refining Estimates of Population Vulnerability

Our assessment of vulnerability does not consider a component of adaptive capacity that is likely to be important for desert bighorn. The ability of individuals of a given genotype to modify their phenotype in response to environmental conditions, known as phenotypic plasticity, is an important mechanism for climate adaptation that can vary among individuals and populations, but it is not well understood and can be difficult to measure (Merilä and Hendry, 2014; Nicotra et al., 2015). Of the in situ mechanisms for adaptation to climate change, plastic responses may be at least as important as evolutionary responses for many species (Hoffmann and Sgro, 2011), although these two mechanisms can be difficult to distinguish (Merilä and Hendry, 2014). Experimental or observational studies of plasticity could provide useful information for prioritizing vulnerable desert bighorn populations (e.g., Renaud et al., 2019).

We avoided using correlative models in our assessment of climate vulnerability because they have several limitations with respect to desert bighorn sheep. Correlative models do not realistically reflect dispersal limitations that influence the geographic distributions of bighorn populations. They may be difficult to apply at the subspecies level because occurrence may be linked with different (and unknown) climatic variables in different parts of the highly variable historic range. Correlative models also typically have not accounted for biotic interactions or changes in vegetation type, which could be critically important for bighorn sheep given that predicted shifts from forest to shrub communities in the southwestern U.S. (Williams et al., 2010) could create new habitat for bighorn. Recent research, however, shows promise of overcoming some of these limitations. For instance, a study of forest bats by Razgour et al. (2019) incorporated information on local climate adaptations from genotype-environment association analysis directly into ecological niche models, and then combined this with a landscape genetic analysis to infer how dispersal barriers are likely to influence potential for evolutionary adaption to climate change. Innovative approaches such as that hold great promise for producing more useful and comprehensive climate vulnerability assessments. In the meantime, less complex approaches such as ours may serve wildlife managers constrained by the amount of time and resources allocated to researching climate vulnerability for any particular species or population.

Management Implications

Maintaining genetic diversity of, and gene flow among, desert bighorn populations will be critical for facilitating evolutionary adaptation to climate change, and actions to preserve or restore connectivity could be effective tools for achieving this objective. A previous simulation study of gene flow for desert bighorn sheep in parts of the study area (Creech et al., 2017) suggested that the spread of adaptive genetic variation is strongly influenced by habitat configuration, with faster spread occurring in areas less fragmented by natural and anthropogenic barriers to dispersal. Reductions in connectivity and gene flow due to anthropogenic barriers such as interstate highways can sometimes be addressed with wildlife crossing structures, although monetary costs can be very high (Corlatti et al., 2009). For example, overpasses for desert bighorn sheep constructed along Highway 93 in the Black Mountains of Arizona facilitated > 1,700 crossings in the first three years after construction (Gagnon et al., 2013). Crossing structures in other areas, such as along Interstates 10, 15, and 40 in southern California, could greatly enhance regional connectivity (Creech et al., 2014), and there may also be opportunities to improve connectivity by modifying existing infrastructure (e.g., highway fencing around underpasses). In some cases, bighorn sheep may eventually discover routes over or under such barriers on their own (Epps et al., 2018). Where connectivity challenges cannot be addressed through barrier mitigation actions, periodic translocation of individuals from nearby populations could be used to provide gene flow (Weeks et al., 2011; Hedrick and Wehausen, 2014).

Ultimately, any management action that promotes large population sizes should help to reduce vulnerability by preserving genetic diversity (Hedrick, 2011). This could include actions that target stressors, whether climatic or non-climatic, such as maintaining or adding artificial water sources (Dolan, 2006; Longshore et al., 2009), controlling predator populations (Wehausen, 1996; Ernest et al., 2002; Rominger et al., 2004), enhancing forage quality (e.g., via prescribed burning; Holl et al., 2004), reducing effects of recreation (Longshore and Thompson, 2013; Wiedmann and Bleich, 2014), or minimizing disease risk by preventing co-mingling with domestic animals (Wehausen et al., 2011). We note, however, that the feasibility of management actions may vary by jurisdiction (e.g., predator control is typically not employed on NPS lands).

Our analysis assumed that connectivity is beneficial to desert bighorn sheep populations because of its critical role in promoting gene flow and maintaining genetic diversity, but connectivity can also facilitate the spread of disease among populations (Simberloff and Cox, 1987; Hess, 1994, 1996). This threat is particularly acute for bighorn sheep, which are susceptible to diseases introduced by domestic livestock beginning with European settlement in the late 1800s (Wehausen et al., 2011). Disease has led to declines across large portions of bighorn sheep range and continues to affect bighorn populations, even within our study area; for instance, respiratory disease outbreaks have been detected in several populations in the Mojave Desert in recent years (Epps et al., 2018). Thus, the risk of disease transmission should be considered before undertaking management actions to increase dispersal among populations. This trade-off between evolutionary potential and disease risk represents one of the fundamental questions that must be clarified for management in the face of both disease and climate change.

Our results underscore the need to maintain native desert bighorn sheep populations in ecologically intact landscapes to facilitate evolutionary adaptation. Unlike other systems where restoration through translocation has dramatically influenced genetic diversity and structure (Whittaker et al., 2004; Malaney et al., 2015; e.g., Jahner et al., 2019), we found that desert bighorn in regions with few anthropogenic barriers to dispersal where populations have persisted without significant reintroduction or augmentation efforts, such as DEVA and GRCA, tended to exhibit high genetic diversity and low isolation. In contrast, populations in landscapes that are more fragmented (i.e., southern California) or heavily influenced by translocations (i.e., southeastern Utah) tended to be less genetically diverse and more isolated. This pattern is consistent with previous studies that have found low genetic diversity or fitness in reintroduced bighorn sheep populations (Whittaker et al., 2004; Wiedmann and Sargeant, 2014). Fortunately, areas that support some of the most genetically diverse and connected populations in our study area are also predicted to have relatively low climate change exposure (e.g., DEVA and GRCA).

The climate vulnerability results from our analysis should help NPS and other natural resource management agencies make more informed decisions about allocating resources among desert bighorn populations, but they do not establish an optimal management strategy. Should efforts be focused on populations with the greatest climate vulnerability in hopes of preserving populations throughout desert bighorn range? Or should managers dedicate resources primarily to bolstering populations that appear most likely to persist in the face of climate change? Current data and analysis cannot fully inform these decisions, and management priorities will also depend on a variety of factors not considered here. For instance, the contributions of genetically unique populations to the evolutionary potential of the species, the susceptibility of connected populations to disease outbreaks, and the public recreational value of populations that provide hunting and wildlife watching opportunities could all influence bighorn management strategies.

We believe the climate vulnerability assessment approach demonstrated here could be applied to the management of many wildlife species other than bighorn sheep. Species that exist in discrete populations due to naturally patchy habitat distribution or the existence of anthropogenic dispersal barriers are likely better suited to this approach than are species that are continuously distributed throughout their range. For well-studied species, previous research identifying specific environmental and climatic factors that influence exposure, sensitivity, adaptive capacity, or landscape connectivity could be used to further tailor this approach. However, vulnerability assessments for poorly studied species may need to rely more on species-neutral measures of climate vulnerability components (e.g., climate velocity or coarse-filter connectivity models based on landscape naturalness). Increasing use of genetic techniques in wildlife management and conservation (Pierson et al., 2016; Smith et al., 2016) suggests that genetic data needed to apply our approach may already exist for some species. Given the capacity for genetic data to expand our understanding of climate change vulnerability, acquiring these data for additional species and populations should be a priority.

Data Availability Statement

The datasets generated and analyzed for this study can be found in the FigShare respository at: https://doi.org/10.6084/m9.figshare.c.4998887.

Ethics Statement

Ethical review and approval was not required for the animal study because the study used non-invasive genetic sampling methods.

Author Contributions

TC, CE, and RM contributed to conception and design of the study. KL was instrumental in the development of the study that added data for 10 populations in the LAKE region. TC, CE, JW, JJ, KL, BH, WS, and RM performed the field work. RC, CE, JJ, and JW performed the lab work. TC wrote the first draft of the manuscript. TC, CE, and RC wrote sections of the manuscript. All authors contributed to manuscript revision and read and approved the submitted version.

Funding

Funding for this research was provided by a Provost’s Distinguished Graduate Fellowship from Oregon State University and NPS’s Climate Change Response Program (PMIS 162673). Sampling and data development in southern Nevada, adjacent Arizona, and southeastern California was funded by the Bureau of Land Management (Agreements FAA010017 and FAA060067) and NPS through LAKE (Coop Agreement H8R07060001, Task Agreements J8R07070011 and J8R07070012). Funding to NPS, LAKE was from the Southern Nevada Public Lands Management Act. The first sampling and data development in the DEVA region in California from the White and Inyo Mountains to Death Valley was funded by the California Department of Fish and Wildlife. Funding for sampling in GRCA was provided by NPS’s Natural Resources Preservation Program (PMIS 157851) and the Grand Canyon Conservancy.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Acknowledgments

We thank S. Ratchford, S. Holt, L. Manning, C. Crow, G. Lenhart, M. Miller, J. Spaak, M. Lamb, and L. Pilkington for assistance with field logistics and sample collection. C. Peterson and M. Corbin helped with genotyping samples and J. Greenwood assisted with extractions. M. Buchalski provided advice on genetic structure analyses.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2020.00279/full#supplementary-material

References