The Silver Bell Mountains of southern Arizona are home to an isolated population of Mexican desert bighorn sheep (Ovis canadensis mexicana) (Jansen et al., 2006). There are few nearby populations of bighorn sheep in south-central Arizona, the nearest is more than 50 km away, separated by inhospitable habitat between suitable mountain ranges, and no immigration has been recorded into the Silver Bells. The Silver Bell Mountains population is the last remnant of an endemic bighorn population complex that once included the Santa Catalina, Santa Rita, and Rincon Mountains in southeastern Arizona (Arizona Game and Fish Department, pers. comm., 2017). Though bighorn sheep were historically ubiquitous across much of the mountainous West, they suffered catastrophic declines in the late 1800s and early 1900s (Valdez & Krausman, 1999). A combination of unregulated harvest, the introduction of livestock-related diseases, and habitat alteration post-European colonization drove many bighorn sheep populations into extirpation and left fragmented remnants of formerly vibrant metapopulations in its wake (Valdez & Krausman, 1999). Though the story may sound familiar, at first glance, the population in the Silver Bells presents an interesting twist on this time-old tale: this isolated population is rapidly expanding after decades of remaining consistently small.
As with all bighorn sheep populations in the state, the Arizona Game and Fish Department (AZGFD) has closely monitored the demographics and health of this population for decades. Until recently, the number of individuals in the Silver Bell bighorn sheep population had not fluctuated dramatically, at least since 1981 when monitoring began (Fig. 1). Annual aerial surveys showed a relatively stable population from 1981 to 2003, with the number of individuals observed fluctuating between 24 to 50 sheep with few exceptions (Fig. 1). In December 2003, outbreaks of infectious keratoconjunctivitis and contagious ecthyma reduced the population approximately 25% through 2008 (Jansen et al., 2007). Subsequently, the population began to slowly increase before rapidly expanding after 2011 (Fig. 1). The number of individuals observed during surveys increased from 34 sheep in 2011 to 133 and 168 in 2014 and 2015, respectively (Fig. 1). Surveys were not conducted in 2012 and 2013, due to changes in AZGFD’s statewide sheep surveying protocols. Although some previous estimates of biological potential for growth in bighorn sheep suggest a doubling of the population in four to five years (Gross, Singer & Moses, 2000; McCarty & Miller, 1998), the number of sheep observed in the Silver Bell Mountains population indicated a more dramatic increase in a four-year span.
Based on its history as an endemic population that underwent a small, albeit recent, disease-induced bottleneck in 2003, we hypothesize that a rare migration event into the population could explain the rapid expansion seen in the Silver Bell Mountains by both increasing abundance and by boosting reproductive vigor through a genetic rescue (Tallmon, Luikart & Waples, 2004; Hedrick, Adams & Vucetich, 2011; Whiteley et al., 2015).
In general, small populations are inherently more vulnerable to change in the form of environmental catastrophes, demographic changes, or stochastic environmental changes, than are large populations (Lande, 1988; Sæther et al., 2005; Willi, Van Buskirk & Hoffman, 2006). Whether isolated naturally or due to anthropogenic forces, these small populations are also more vulnerable to genetic drift, resulting in the loss of genetic diversity, fixation of alleles within populations, and reduced evolutionary potential (Frankham, 1996; Kimura, 1955; Otto & Whitlock, 1997; Wright, 1931). Additionally, inbreeding, which lowers genetic diversity by increasing homozygosity, is inevitable when the population becomes sufficiently small and isolated (Charlesworth & Charlesworth, 1999).
Small populations can often benefit in a number of ways from immigration, the introduction of individuals from other populations (Tallmon, Luikart & Waples, 2004; Whiteley et al., 2015). Immigration boosts the number of individuals, making the population more robust while also having positive genetic consequences. Specifically, genetic rescue is an increase in population fitness, as measured through demographic growth, caused by the introduction of novel alleles through immigration (Whiteley et al., 2015). Migrants may introduce novel alleles which can mask the negative effects of deleterious alleles that have built up over time (Tallmon, Luikart & Waples, 2004). The offspring of these migrants may see a boost in fitness related to higher heterozygosity, known as heterozygote advantage (Tallmon, Luikart & Waples, 2004). There is evidence that the arrival of even a single immigrant can be enough to cause the rapid spread of new alleles and exponential population growth, halting the effects of loss of heterozygosity and inbreeding depression (Vila et al., 2003). In the past, genetic rescue has improved longevity and fitness in bighorn sheep populations based on empirical evidence (Miller et al., 2012).
Despite the isolation, bighorn sheep have been known to disperse into neighboring habitats (Singer et al., 2000), making it possible that immigrant sheep could enter the Silver Bell Mountains population. Though rams and ewes have been shown to engage in natal dispersal with similar frequencies in long, linear mountain ranges (Buchalski et al., 2015), in the sky islands habitat of southern Arizona, females are believed to move less frequently between ranges, making male-biased dispersal more common (Buchalski et al., 2015; Epps et al., 2005; Festa-Bianchet, 1991). In O. c. nelsoni, the estimated maximum dispersal distances for rams and ewes in the Mojave Desert were ten times larger across sloped terrain (164 km and 100 km for rams and ewes respectively) than across flat terrain (16.4 km and 10 km respectively) (Epps et al., 2007; Creech et al., 2014). In other studies, dispersal over more than 60 km has been deemed unlikely (Buchalski et al., 2016). With this in mind, it is possible that the rapid growth of the Silver Bell Mountains population could be the result of immigration from other populations in Arizona.
Thus, we applied a multilocus genetic approach to determine whether immigration can at least partially explain the recent, rapid population expansion observed in the Silver Bell Mountains. If immigration was a contributing factor in expansion, we would expect to see a change in allele frequencies, and likely new alleles, in the Silver Bell Mountains population post-expansion.
We used mitochondrial DNA sequence and nuclear microsatellite markers for genetic analyses of the Silver Bell Mountains bighorn sheep before (2003) and during (2015) the population expansion. We also analyzed a small number of available blood samples from the Gila Mountains (a southwestern Arizona population) and the Morenci Mine (Rocky Mountain bighorn sheep O. c. canadensis) to detect the source of putative immigrants and, more importantly, to serve as comparisons for genetic diversity metrics. The Silver Bell Mountains population is of management concern for AZGFD (AZGFD, pers. comm., 2017), and we thus provide the only published population genetic analyses specifically focused on this population of concern.
Materials and Methods
Sampling and DNA extraction
We obtained 62 bighorn sheep blood samples collected by AZGFD (Fig. 2). Of these, 52 were from the Silver Bell Mountains population (35 from the 2003 capture, when AZGFD was treating bighorn sheep for disease, and 17 from a 2015 capture, after the population had crashed and subsequently increased), five from the Gila Mountain Range, collected in 2016 from south of Highway 8 in southwestern Arizona, and 5 Rocky Mountain bighorn sheep samples collected from Morenci Mine in 2014. We used the FlexiGene DNA blood extraction kit (Qiagen, Inc., Valencia, CA, USA), and followed the manufacturer’s protocol to extract DNA from whole blood stored in EDTA tubes.
We used the primers L15712 for the forward sequence and BETH for the reverse sequence (Epps et al., 2005) to amplify a 515 base pair fragment of the mitochondrial control region. We ran 20 µL PCR reactions with the following conditions: 1×PCR Buffer (Invitrogen™, Thermo Fisher Scientific Inc. Waltham, MA, USA), 0.16 mM dNTPs, 0.05% bovine serum albumin (Sigma-Aldrich, St. Louis, MO, USA), 1.9 mM MgCl2, 0.4 uM of each primer, 0.8 units of Taq DNA polymerase (Invitrogen™), and 10 ng of extracted DNA. Amplification conditions were as follows: 94 °C for 7.5 min, 35 cycles of 94 °C for 60 s, 61 °C for 70 s, and 72 °C for 90 s, and a final extension at 72 °C for 7 min. We visualized PCR products by gel electrophoresis on 1.5% agarose gels stained with Ethidium Bromide. We used EXOsap-IT (USB-Affymetrix, Inc.) to purify PCR products prior to sequencing with the following PCR conditions: 37 °C for 15 min followed by 80 °C for 15 min then held at 4 °C. We sent final purified PCR products to the University of Arizona Genetics Core (UAGC) for sequencing in both the forward and reverse directions on an ABI3730 DNA Analyzer (Applied Biosystems, Thermo Fisher Scientific Inc., Waltham, MA, USA).
We used Sequencher v.5.0.1. (GeneCodes, Ann Arbor, MI) to assemble the resulting forward and reverse sequences manually. We visualized chromatograms, aligned and trimmed to a reference sequence (Genbank ID: KP688366) to resolve ambiguous reads, and exported contiguous (FASTA) sequences for analyses. We imported the sequences into MEGA v.7 (Kumar, Stecher & Tamura, 2016) to obtain a final alignment. A 75 base pair repetitive sequence was discovered to be duplicated in one haplotype, RM6. Following Buchalski et al. (2016), samples were reduced to a single copy of the repeat. We trimmed all sequences to 515 base pairs.
We used 15 microsatellite markers, which were split into two multiplexes of 8 and 7 microsatellites each (Multiplexes 1 and 2 from Buchalski et al. (2015)). For 20 uL reactions, concentrations were 1×Invitrogen reaction buffer, 2.5 mM MgCl2, 0.4 mM dNTP’s, 0.12% BSA, and 1 unit of Taq DNA Polymerase (Qiagen Inc., Valencia, CA, USA; see Buchalski et al., 2015 for primer concentrations). We used a touchdown approach for amplification conditions: 94 °C for 15 min, followed by 5 cycles of 94 °C for 30 s, annealing at 63 °C to 59 °C decreased by 1.0° per cycle for 90 s, 72 °C for 1 min, and then 30 cycles of 94 °C for 30 s, 59 °C for 90 s, 72 °C for 1 min, and a final extension at 60 °C for 30 min. We sent PCR products to the UAGC for fragment size analysis in an ABI3730 DNA analyzer (Applied Biosystems TM, Thermo Fisher Scientific Inc., Waltham, MA, USA). We used Genemarker v.2.6.0 (SoftGenetics, State College, PA, USA) to score allele sizes. Three to four samples were included as duplicates on each plate we genotyped to ensure that our allele calls were consistent.
For mtDNA, we exported a final alignment from MEGA v.7 into DnaSP v.5.10.1 (Librado & Rozas, 2009). We used DnaSP to estimate the number of haplotypes, haplotype diversity, and nucleotide diversity. We used the algorithm BLAST (Altschul et al., 1990) to match mitochondrial haplotypes with haplotypes previously recorded in GenBank (Benson et al., 2018). We used PopArt to create a median neighbor joining haplotype network (Leigh & Bryant, 2015).
Three of our microsatellite markers failed to consistently amplify and as a result our final data set was reduced to 12 markers. For the nuclear microsatellites, we used GenAlEx v.6.501 (Peakall & Smouse, 2012) to determine if loci were in Hardy-Weinberg equilibrium, measure allele frequencies, and heterozygosity. Additionally, to adjust for our large disparity in sample sizes between different populations, we used HP-Rare to calculate allelic richness using rarefaction (Kalinowski, 2005). We calculated pairwise- FST between the four different populations using Arlequin 188.8.131.52 (Excoffier & Lischer, 2010). Calculations were based on Weir & Cockerham (1984) and statistical significance was determined using 10,000 permutations and a p-value of 0.05. We used GeneClass2.0 to detect first-generation migrants among each population (Piry et al., 2004). We applied the Paetkau et al. (1995) frequency-based criterion for likelihood computations, a frequency of 0.01 for missing alleles, the Paetkau et al. (2004) resampling algorithm, 10,000 resampled individuals, and a threshold significance of p < 0.01.
Additional analyses were run specifically in the Silver Bell populations. Probability of Identity (P(ID)) (Paetkau et al., 1998) and Probability of Identity assuming siblings (P(ID)sib) (Waits, Luikart & Taberlet, 2001) were calculated and we looked for samples with identical genotypes between time points using GenAlEx (Peakall & Smouse, 2012). Using the R package Rhh (Alho, Välimäki & Merilä, 2010), we calculated homozygosity by locus (HL) and internal relatedness (IR) for each individual (Aparicio, Ortego & Cordero, 2006). Welch’s t-test was used to determine whether or not differences between HL and IR for the two temporal samples were significantly different. Finally, we used ML-Relate to calculate the pairwise relatedness (r) between each individual in each temporal sample (Kalinowski, Wagner & Taper, 2006).
We used the Bayesian clustering procedure implemented in STRUCTURE v184.108.40.206 (Hubisz et al., 2009) to genetically assign individuals to clusters, assuming admixture and correlated allele frequencies between populations. We considered the number of populations (K) from 1 to 10. We ran STRUCTURE with 1,000,000 Markov chain Monte Carlo iterations and 20 independent simulations per K with the first 200,000 iterations eliminated as burn-in. We used the ΔK statistic of Evanno, Regnaut & Goudet (2005) implemented in the CLUMPAK server (clumpak.tau.ac.il; Kopelman et al., 2015) to determine the most appropriate number of genetic clusters.
From the mitochondrial DNA, we observed six haplotypes among the collected samples, all of which were matched to known haplotypes from GenBank (Fig. 3). We found three haplotypes within the Gila population, and two haplotypes from the Morenci Mine population. The remaining haplotype was shared in all samples from the Silver Bell Mountains, from both 2003 and 2015. There were fewer base pair differences between the Silver Bell Mountains population and the Gila population, both of which are O. c. mexicana, than the Morenci Mine population, which is O. c. canadensis (Fig. 3). On average, haplotype diversity was 0.299, and the nucleotide diversity was 0.007 (Table 1).
At the microsatellite loci, the sheep from the Silver Bell Mountains had lower genetic diversity (Table 1). Furthermore, sheep from the Silver Bell Mountains had lower allelic richness (1.54 in 2003 and 1.57 in 2015) than those from the Gila Mountains (2.33) and Rocky Mountain bighorn sheep from Morenci Mine (2.92). The Silver Bell Mountains population showed a much lower unbiased expected heterozygosity (0.19 for 2003 and 0.21 for 2015) than the Gila and Morenci Mine populations (0.43 and 0.52 respectively). FIS values from all populations were near zero, with −0.052 for Morenci Mine, −0.143 for Gila, 0.001 for the Silver Bell Mountains in 2015, and 0.014 for the Silver Bell Mountains in 2003. Additionally, all pairwise FST estimates were significant, except for the comparison between the 2003 and 2015 Silver Bell Mountains sheep (Table 2). We did not detect any statistically significant first-generation migrants with GeneClass2.0.
The two temporally separate sample sets from the Silver Bell Mountains population were similar genotypically, but not identical. The 2003 sheep had two alleles, both at locus MAF214, which did not appear in the 2015 population. Conversely, we did detect one allele in an individual in the 2015 population that was not detected in any of the 2003 sheep sampled. Average pairwise genetic relatedness (r) between individuals was comparable for each temporal sample (0.198 for 2003 and 0.176 for 2015). Additionally, average HL for each temporal sample was measured at 0.56 in 2003 and 0.54 in 2015; the Welch’s t-test suggested that this difference was not significant (p = 0.815). Likewise, average IR for each temporal sample equaled 0.056 in 2003 and then −0.005 in 2015, yet again this was a non-significant difference (p = 0.624). The STRUCTURE plot failed to detect a significant shift in allele frequencies between the 2003 and 2015 Silver Bell Mountains bighorn sheep population as well (Fig. 4). Finally, P(ID) was calculated as 7.6 ×10−3 and 5.1 ×10−3, for 2003 and 2015. P(ID)sib was calculated as 9.2 ×10−2 and 7.2 ×10−2 respectively. Six samples in 2003 had a matching genotype to a sample in 2015, and there were 3 pairs of samples in 2003 that had the same genotype at every locus. It is possible that some sheep were resampled in 2015, but our detection thresholds are too high to rule out the possibility that these were different individuals, especially since we found some identical genotypes in just the 2003 samples.
Our original hypothesis to explain the population expansion seen in the Silver Bell Mountains was that immigration into the population resulted in the observed expansion through both numerically increasing the number of sheep and affecting a genetic rescue. However, we are unable to support this hypothesis. Based on the genetic data, we failed to support findings of substantial gene flow into the population. We found no new mitochondrial haplotypes in the 2015 sheep from the Silver Bell Mountains and only a single new microsatellite allele, found in one individual. While the lack of novel mitochondrial diversity is less surprising, based on the general pattern of male-biased dispersal in desert bighorn sheep, the microsatellites would pick up novel alleles introduced by male immigrants. It is likely that the one new microsatellite allele was simply not measured in the 2003 population due to its low frequency and the sample sizes involved. Interestingly, this rare allele was prevalent in the Rocky Mountain bighorn sheep from the Morenci Mine site; yet no other evidence from the mitochondrial or nuclear genetics supports immigration from the Morenci Mine population. GeneClass2.0 did not highlight a single sample in the 2015 population as a likely first-generation migrant, though this analysis is limited by only having the two potential source populations for migrants.
The Silver Bell Mountains population of bighorn sheep did show some genetic signals of being small and isolated. Genetic diversity was much lower in the Silver Bell Mountains than seen in either the Gila population or the Morenci Mine population; there was more diversity present in just those five samples than in the entire Silver Bell Mountains. This finding, in spite of our low sample size (only five samples in each of those two populations) is worrying for the Silver Bell Mountain sheep population. However, concerns about the Silver Bell Mountains population being potentially inbred do not seem warranted based on our data. FIS, HL, IR and r all provide some support that this population, despite its small size for many generations seemingly without immigration, has managed to mostly avoid inbreeding. For example, other isolated, bottlenecked populations of desert bighorn have much higher estimated inbreeding coefficients: Red Rock, AZ at 0.390 and Tiburon Island, MX at 0.288 (Hedrick, 2014) compared to the Silver Bell Mountains at 0.001. While we did not detect genetic signals of inbreeding, immigration would still likely be beneficial to the genetic health of the Silver Bell Mountains population, especially in the face of its low genetic diversity and continued isolation.
While we did not detect gene flow into the Silver Bell Mountains, we cannot definitively reject the possibility that immigration has occurred, especially in recent years. The lack of new mitochondrial haplotypes only rules out immigration from females with a new haplotype, meaning males, or females with this same haplotype, could have entered the population. The haplotype observed in the Silver Bell Mountains is not unique and it has been observed elsewhere in Arizona and California (Buchalski et al., 2016; Epps et al., 2005). The lack of novel microsatellite alleles, on the other hand, is strong evidence for the lack of gene flow into this population. However, if gene flow into the Silver Bell Mountains population occurred from a population with a very similar genetic makeup it might be undetectable with our methods but may be differentiated enough to provide a slight increase in diversity and improved fitness. Additionally, functional variation could have been introduced into the Silver Bell population and played a part in the recovery. We would not have been able to detect functional variation in the nuclear genome, as we only used neutral microsatellites. An increase in MHC diversity could have played a part in the expansion, as this population had a recent history of disease-induced bottlenecks.
Our inability to detect novel alleles or changes in allele frequencies in the Silver Bell Mountains population before and after the population boom could alternatively be a function of timing. Introgression into the population, especially if from a closely related population, would see the novel alleles slowly filter through the populations (Johnson et al., 2010). By sampling in 2015, the novel variation could simply be too low in frequency to have been detected with our sample size and number of markers. For example, when Florida panthers (Puma concolor) were rescued by an introduction of pumas from Texas, it took a minimum of 3 years, about one full generation (Logan & Sweanor, 2001), before Texas alleles became detectable in the Florida population (Johnson et al., 2010). With that being said, generation time for bighorn sheep is usually considered to be about 6 years (Coltman et al., 2003; Johnson et al., 2011). Generally, once a ewe reaches the age of two or three, they birth one or, rarely, two lambs per year (Monson & Sumner, 1980). Rams, on the other hand, are generally not successful breeders until the age of four (Monson & Sumner, 1980). If we assume that two generations have occurred during that time frame, it is possible that this is too few generations post-expansion to see changes in allele frequencies at only 12 microsatellite markers. However, Epps, Crowhurst & Nickerson (2018) detected the influx of migrants into desert bighorn sheep populations over just two generations (2003–2015) using a similar number (15) of microsatellite markers. They found significant changes in genetic structure due to immigration and were able to determine which sheep were immigrants and the populations from which they originated (Epps, Crowhurst & Nickerson, 2018). In order to detect potential evidence of immigration that we could have missed with our study design, a future study could employ some combination of more microsatellites markers, larger sample size post-expansion, and sampling generations into the future. A study of adaptive variation, whether that be through genome wide SNPs or specifically the MHC, could elucidate potential differences in functional genes that have fitness consequences.
Alternatively, if no immigrants entered into the population, then the Silver Bell Mountains bighorn sheep population could have expanded due to a change in some environmental variables. Precipitation (Bender & Weisenberger, 2005), escape cover (Dunn, 1996), predation (Hayes et al., 2000), forage quantity and quality (DeYoung et al., 2000), and disease (Gross, Singer & Moses, 2000) all have been shown to have effects on growth of bighorn sheep populations, and thus some potentially undetected change to these environmental conditions could have, in theory, driven the growth instead. At full biotic capacity, bighorn sheep are believed to be capable of doubling their population size every 4–5 years (McCarty & Miller, 1998), yet based on standardized, repeated aerial helicopter surveys, the number of sheep observed in the Silver Bell Mountains nearly quadrupled in that same time frame. The observed population growth is exceptional, even accounting for wide margins of error due to survey methodology.
Population genetic comparison of two time points and two potential source herds using 12 microsatellite loci and a mitochondrial marker did not support immigration as an explanation for population expansion. We hypothesized that we would find an empirical example of unassisted genetic rescue in the wild, yet we failed to find evidence that the molecular benefits of immigration attributed to this dramatic population growth and expansion. Though we did not explain this phenomenon, we do offer the first published genetic characterization of this population of concern. Continued monitoring (abundance, habitat and environmental variables, health characteristics such as disease, etc.) is more important than ever in the face of this intriguing population trend, especially as we demonstrated that the Silver Bell Mountains population has reduced genetic diversity. Determining whether the population is continuing to grow, is at carrying capacity, or if it will disperse to adjacent mountain ranges could be important for determining why the population growth occurred in the first place. Based on our mitochondrial and nuclear genetic analyses, this population expansion does not appear to be directly linked to a genetic or demographic effect of immigration, though increased sampling efforts with more individuals, potential source populations, and non-neutral genetic markers over a longer time period could provide a clearer picture. An environmentally-linked driver such as changes in habitat, predation, or something else entirely has likely resulted in this population expansion, though we cannot emphatically reject our original hypothesis of genetic rescue.
Delta K graph for STRUCTURE
This is the Delta K graph used to select which K values best fit the model. We used the K statistic of Evanno, Regnaut & Goudet (2005) implemented in the CLUMPAK server (clumpak.tau.ac.il; Kopelman et al., 2015) to determine the most appropriate number of genetic clusters.