Grand Canyon is an isolated, 2.4 km-deep (1.6 km average), 28 km-wide (16 km average) canyon on the southern Colorado Plateau in northern Arizona that spans nearly 450 river kilometers and is located at the geographic limits of the Great Basin, Mojave Desert, and Sonoran Desert. Grand Canyon consists of plateaus and other comparatively flat topography incised by shallow drainages on the plateaus and steep canyons leading into the Colorado River gorge. On the rims, natural surface water is uniformly scarce, but sparsely distributed springs and seeps serve as perennial sources of water in the inner canyon. Collectively, the area includes approximately 6,750 km² of relatively steep terrain below the rims of Grand Canyon. Our study area primarily focused on 364 km of the river, which started at Lees Ferry (river mile (RM 0) and ended at Diamond Creek (RM 226) but included light-intensity sampling and observations along another 85 km of river, terminating at Pearce Ferry on Lake Mead (RM 279). We use RM (Gushue, 2019) to denote location (instead of km) due to the widespread, universal use of this nomenclature for the Colorado River in Grand Canyon. Similarly, given the meandering nature of the river, we use the term river left (RL) or river right (RR) to denote river side as when viewing downstream from RM 0 to RM 226.

Hoover and Glen Canyon Dams on the Colorado River create the nation’s two largest reservoirs, Lake Mead and Lake Powell, which roughly form the lower and upper ends, respectively, of Grand Canyon. The Grand Canyon is contained within the park and lands of the Hualapai Indian Tribe. The Havasupai Tribe, the Navajo Nation, and a mosaic of other federal land management entities. Elevations range from 600 meters along the canyon bottom, where the Colorado River flows, to over 2,400 meters on north rim of Grand Canyon ( Figure 1 ).

Grand Canyon is geologically divided into the eastern and western basins (Billingsley and Hampton, 1999), separated by the steep, cliff-bound Muav Gorge that serves as the primary upstream and downstream barrier for dispersal of several plant, invertebrate, and vertebrate taxa (Miller et al., 1982; Phillips et al., 1987; Stevens and Polhemus, 2008). River flow rates change depending on daily discharge, which varies in response to hydropower demand and river reach, averaging 226–340 cubic meters per second (8,000–12,000 cubic feet per second). The rims of Grand Canyon routinely experience snowfall in winter. In contrast, the lower levels of Grand Canyon near the river rarely drop below freezing (January minimum of 2.4°C), are hot during the summer (mean July maximum of 42.0°C), and receive an average of only 225 mm of precipitation. Precipitation occurs throughout the year except during the strong spring drought during May and June, with highest precipitation during a summer monsoonal period from July through September.

Below the canyon rim, Mojave and Sonoran Desert scrub communities dominate the lowest elevations, rising to more cold-tolerant species indicative of Great Basin Desert communities. Along the riparian zone, vegetation structure has been modified since the completion of the Glen Canyon Dam, which generally has prevented floods that used to scour the river channel and limit plant colonization (Sankey et al., 2015). Stabilized water flows have since resulted in dense vegetation along the shoreline (Sankey et al., 2015), frequently dominated by nonnative tamarisk (Tamarix ramosissima). Above the shoreline zone, vegetation composition varies with aspect and solar radiation (Stevens, 2012), but is dominated by species including honey mesquite (Prosopis glandulosa var. torreyana), catclaw acacia (Seneglia greggii), brittlebush (Encelia farinosa), creosote (Larrea tridentata), and numerous yuccas and cactuses. More than 300 plant species have been described within Grand Canyon.

Bighorn in Grand Canyon are most commonly observed along the river corridor. Since the 1950’s, anecdotal observations of bighorn along the river were periodically collected by park and commercial (tourism) river companies; for decades, this was the primary source of information about the population. In 2000, park staff began to work with private and commercial river guides to report bighorn sightings in a standardized way, reporting sex/age class, group composition, river mile, and river side. We generated location coordinates at the river mile intersection with the modeled shoreline for 20,000 cfs (566 cms) flows. We used 2000–2018 bighorn demographic data collected by park, commercial, and private river trips across all seasons to inform numbers and distribution of bighorn along the Colorado River through the park. Data collected between 2000–2010 provided a more comprehensive assessment of bighorn distribution along the river than was previously available and helped guide our study beginning in 2010. Bighorn congregate along the river during the breeding season (August–October) and are most concentrated during that time. Thus, in 2010, we began counting and classifying bighorn during September bighorn-specific river trips. Three to four experienced observers surveyed for bighorns on each trip, but depending on the timing of the observation relative to the position of the raft and the ability of the boat operator to stop in that location, the time for counting and classifying a group of bighorn varied, and observers worked together to document all sightings. Detectability of bighorns from the river via visual observations may vary based on time of day, temperature, terrain ruggedness, and activity level, and robust population estimates from this approach would require independent double observations or large numbers of marked animals. However, this can provide insight to distribution over time and classifications provide an index to reproduction with enough observations. In 2011, we broadened the study to include information on disease surveillance and genetic structure. River trips usually consisted of four biologists and one boat operator in a single 22-foot motor raft that allowed us to hold position except in the largest rapids, but did not allow efficient upriver movement. In 2011 we had 2 boats. In total, we made seven annual river trips between 2010 and 2016, with each lasting approximately two weeks. Demographic counts continued through 2018, by which time all collared bighorn either dropped their collars or died.

We recorded all observations of bighorn with symptoms of respiratory disease and mortality. We classified bighorn as sick if they expressed clinical signs of respiratory complications including nasal discharge, coughing, and labored breathing, and we investigated sick or dead bighorn reported by commercial river trips when feasible from 2010–2018 ( Supplementary Tables A1, A2 ). To determine current and prior exposure to disease pathogens, we collected antemortem diagnostic samples from bighorn captures from 2013–2016 (see below), and postmortem diagnostic samples were collected from recent mortalities when feasible from 2011–2018. For captured bighorn, approximately 20 mL of blood was collected into serum-separating tubes as well as a 5 ml tube with EDTA for whole blood preservation for DNA work and kept on ice while in the field; after finishing each river trip, serum-separating tubes were immediately spun to extract serum which was then frozen. For captured bighorn and relatively intact carcasses, nasopharyngeal swabs were used to collect nasal and tonsillar samples, stored in tubes with tryptic soy broth (TSB) buffer in glycerol, and frozen. Serum and swabs were analyzed by Washington Animal Disease Diagnostics Laboratory (WADDL) at Washington State University: swabs were tested by PCR for evidence of infection by Mycoplasma ovipneumoniae, a key respiratory pathogen, while ELISA was used on serum samples to test for presence of antibodies to M. ovipneumoniae and was used as a metric for exposure. Strain typing for M. ovipneumoniae was conducted on a subset of samples by WADDL as described by Kamath et al. (2019). On a subset of samples, the Colorado State Veterinary Laboratory conducted PCR tests for other respiratory pathogens, including Biberstenia trehalosi and Mannheimia (Pasteurella) haemolytica. Although respiratory pneumonia in bighorn is a polymicrobial disease, M. ovipneumoniae has been identified as a primary causal agent (Cassirer et al., 2018).

We fitted GPS collars to bighorn to characterize habitat use and movement and to inform the population estimation model. We captured bighorn from 2010–2016 along the Colorado River in Grand Canyon using river- and ground-based darting techniques during the same river trips described above. River-based captures involved darting bighorn from the boat in eddies or slow-moving areas of the river where distance to the animal was less than 30m. Ground-based captures involved approaching and darting bighorn from the shoreline or farther from the river if darting could not be accomplished safely from the boat. Capture efforts were concentrated during fall when bighorn tend to congregate along the river, and to reduce temperature- and gestation-related stress during captures. We attempted to capture animals equally on either side of the river, and spatially distributed throughout the river corridor, except where suitable habitat or shoreline access was absent.

Bighorn were immobilized with BAM, a pre-mixed combination of Butorphanol (B), Azaperone (A), and Medetomidine (M), with a dosage of 0.50 mg/kg B, 0.17 mg/kg A, and 0.20 mg/kg of M. Bighorn were aged and classified based on annual horn rings and time of year (Elbroch, 2006). We fitted bighorn with GPS-satellite uplink (Telonics TGW-4583) or GPS-direct download (Telemetry Solutions Quantum 5000) collars programmed to collect GPS fixes every four hours and equipped with a mortality sensor and a programmable drop-off mechanism. Locations of marked bighorn were transmitted via Argos or Iridium satellite uplink or downloaded directly via UHF antenna during aerial surveys. Bighorn captures were approved by the National Park Service (NPS) Institutional Animal Care and Use Committee (IMR_GRCA_Holton_Bighorn Sheep) and the park (GRCA-2011-SCI-0018, GRCA-2012-SCI-0020, GRCA-2013-SCI-0057, CE-2010-2015 PEPC-21840, CE-2016-2021 PEPC-62989).

We created two genetic datasets to address different purposes. For the primary dataset, used for population estimation, we obtained DNA via non-invasive sampling of fecal pellets during annual river trips in September–early October from 2011–2015. Where accessibility allowed, we collected fecal pellets by 1) searching at locations where bighorn routinely access water, such as the mouths of side drainages, 2) searching where bighorn were observed from the boat, 3) searching where bighorn tracks were observed at the river’s edge, and 4) opportunistic searching at river camps. In the first 3 years, searches of shorelines and adjacent areas (canyons and hills) for fecal pellets was a mix of systematic stops at any promising and accessible location or in response to observations of bighorn or tracks, but we also identified 74 high-priority planned search areas after 2013. These sites were distributed more evenly along the canyon than initial sampling and were either at side canyon–river junctions likely to funnel bighorn or at places with consistent observations in initial sampling years. We stopped at 43 of these in 2014 and 56 in 2015. Samples ranged in age from fresh to those estimated to be up to 2 weeks old based on pellet condition and tracks. Pellet condition, date, UTM location, and river side were noted for each sample, which were collected by gloved hands or using sticks to prevent contamination and placed into a paper envelope. Samples collected wet from freshness or light rain were dried by placing envelopes in direct sun, or by placing envelopes into a pot of dry sand and heating on a camp stove until no moisture was detectable through the envelope; samples were then stored in the laboratory in dry, dark conditions at room temperature.

We used the second genetic dataset for estimation of genetic structure, gene flow, and river crossings (hereafter, “combined” dataset). This dataset was created by merging the primary dataset for population estimation (9 microsatellite loci for 1,250 samples; see Results) with data from samples collected from 2011–2013 and genotyped at 16 microsatellite loci including the 9 used for population estimation (see details in Creech et al., 2017, Creech et al., 2020). Like the primary dataset, that effort relied on non-invasive sampling of fecal pellets, but also included a small number of tissue and blood samples from live captures, euthanized animals, or carcasses found in the field.

We used a modified AquaGenomic Stool and Soil protocol (MultiTarget Pharmaceuticals LLC, Salt Lake City, UT) to extract DNA from material scraped from the surface of fecal pellets. The primary dataset (population estimation) was genotyped at a single panel of 9 microsatellite loci and a sex-identification marker to allow individual identification (see Pfeiler et al., 2020 for details on markers and amplification conditions). The primer pair used for sex identification amplifies the amelogenin gene located on both the X and Y chromosomes (Yamamoto et al., 2002); in bighorn, the Y chromosomal fragment exhibits a 44-base-pair (bp) deletion relative to the X chromosome (214 bp vs. 258 bp). The Creech et al. (2017) dataset was genotyped at 16 dinucleotide microsatellite loci in three multiplex PCRs of 4–6 loci using a Qiagen Multiplex PCR kit (Qiagen, Valencia, CA) and included numerous samples also used in the primary dataset. We submitted PCR products to the Center for Quantitative Life Sciences at Oregon State University for analysis on an ABI 3730 capillary sequencer (Applied Biosystems [ABI], Foster City, CA, USA) and used GENEMAPPER (version 4.1; ABI) to score genotypes. Each sample was amplified in at least three replicate PCRs to generate consensus genotypes. Samples that produced ≥3 alleles at any locus were eliminated. Samples that produced partial genotypes at ≥50% of the microsatellite loci were subjected to 3 more attempted amplifications. For a genotype to be accepted, each allele in a heterozygous genotype had to be observed twice, while the single allele in a homozygous genotype had to be observed 3 times.

We used CERVUS version 3.0.3 (Kalinowski et al., 2007) to identify duplicate genotypes, i.e., those derived from different fecal samples from the same individual bighorn, and used GIMLET version 1.3.3 (Valière, 2002) to estimate genotyping error rates (false allele occurrence rate and allelic dropout rate). We used a maximum probability of identity (P_ID) of 0.01 for unrelated individuals and 0.05 for full siblings (P_IDsibs), and on that basis determined the minimum number of loci (working from least to most informative) needed to be amplified for a genotype to be used. Duplicate genotypes in the primary dataset (9 loci) were identified for population estimation analysis. For the combined dataset, we used genotypes of unique individuals from the primary and Creech et al. (2017) datasets, checked for duplicates after combining them, and removed all but one sample from any individual detected multiple times, retaining the most complete genotypes or selecting at random. Any individual detected on both sides of the Colorado River at either stage was noted as evidence of river crossings. Thus, the combined dataset for evaluation of genetic structure consisted of unique individuals genotyped at either 9 or 16 loci. As Creech et al. (2017) tested these markers and found no evidence of linkage disequilibrium in this system after controlling for genetic structure, we used GENEPOP version 4.2 (Raymond and Rousset, 1995) to test for deviations from Hardy–Weinberg proportions on the combined dataset for genetic structure, after assigning individuals to populations, but did not reevaluate linkage disequilibrium (see below).

We used the combined dataset for analyses of genetic structure and gene flow. Because the distribution of bighorn within the park is discontinuous, but clear geographic boundaries other than the Colorado River were largely not apparent, we used STRUCTURE (Pritchard et al., 2000) with no prior information on population or sampling location to describe clusters of genetically similar individuals. We used a burn-in of 50,000 iterations followed by 100,000 sampling iterations on 5 replicate runs, reviewing consistency across runs to determine that these lengths were adequate. We determined the largest number of clusters diagnosable in the dataset using the “high assignment criterion” (HAC) approach, in which the largest supported value of K is that for which all clusters have at least 1 individual assigned to them at high confidence (Epps et al., 2018). We used CLUMPP to average assignment probabilities (q values) across replicate runs and used a cutoff of q > 0.90 when applying the HAC approach. We also plotted ΔK values (Evanno et al., 2005) against K to evaluate whether peaks at higher K values occurred.

After determining the largest supported value of K, we mapped average assignment probabilities for each sampled individual, choosing at random among multiple locations at which any individual was detected if observed more than once. We then defined populations spatially based on concentrations of samples assigned to coherent clusters, defining them by side (RR or RL) and RM. We reviewed population assignments, observations of bighorn (2011–2018), locations of all genotyped samples, and GPS collar data to guide choices of divisions along each river side, seeking to identify gaps in bighorn distribution where assignment values differed. We used those populations as the basis for subsequent population-based analyses. We tested for migrants among those spatially defined populations using the average q values from STRUCTURE analyses, defining migrants as any animal assigned at q > 0.9 to another cluster.

We evaluated genetic structure among those populations by estimating F _ST using GENEPOP (Raymond and Rousset, 1995). We also tested for evidence of additional structure or other processes such as selection by testing for deviation from Hardy–Weinberg Proportions using the probability test in GENEPOP by population and locus, using default settings for Markov chains.

Finally, we estimated recent gene flow (previous generation) using BIMr (Faubet and Gaggiotti, 2008) among those spatially defined populations. We tested two environmental factors that might influence migration among populations: 1) whether populations were on the same or different sides of the Colorado River and 2) distance, measured along the river from the center of each spatially defined population, where center was defined as the middle point of that stretch of river. We used 10 replicate runs with 20,000 iterations for burn-in and 20,000 iterations for parameter estimation, with other model settings as default. We evaluated posterior model probabilities for 4 models: 1 for each environmental factor, 1 with both factors, and 1 with both factors and an interaction term. We used mean migration rates for each population pair estimated from the run with the lowest Bayesian deviance as a comparison with patterns of migration inferred from STRUCTURE.

Estimating population size posed a particular challenge. Many traditional methods to estimate population size require even sampling across a population and multiple sampling events within a short time window. In the last two decades, spatial capture–recapture methods that allow for a heterogeneous distribution of sampling and integration of data sources have expanded options to estimate population size and understand the role of habitat on distributions simultaneously (Borchers and Efford, 2008; Royle et al., 2013a, Royle et al., 2013c). Ongoing developments now support estimates for wildlife distributed along linear features, such as rivers (Royle et al., 2013b; Fuller et al., 2016; Sutherland et al., 2018), and the use of a single annual sampling session (Morin et al., 2018; Schmidt et al., 2022).

We used a spatial capture–recapture (SCR) model to estimate population size of Grand Canyon bighorn (e.g., Royle and Young, 2008). SCR models detection using a function where the probability of detection declines as a function of distance between the coordinates of a trap, j, and an activity center, s, of an individual animal, i. The activity center is a location, with coordinates, that is analogous to a home range center for a particular period of time, which here likely comprises a few weeks before sampling as the DNA from pellets decays over time. The activity center is a latent variable, that is estimated in this model based on a grid of points (the state space) that describes the potential activity centers considered possible during modeling. We used the half-normal detection function:

This model thus contains 3 parameters, namely density and 2 parameters describing detection, the baseline encounter probability, p ₀, and the spatial scale parameter, σ, which determines the rate of decrease in detection probability with the distance between the trap location, x_j , and the activity center for an individual, s_i . The incorporation of a spatial model for detection allows estimation of detection parameters even with a single detection session per year. We used the approximate place where the boat stopped for the search for pellets to define the trap location and assigned detections of individuals to the trap closest to the pellet collection coordinates. Potential activity centers (the state space) had a spacing of 1 km and included a buffer of 5 km from the river where sampling occurred but extended farther in places to include the entire park.

We summarized sample sizes and calculated the mean maximum distance and the maximum distance among individual detections, confirming >30% of recaptured animals had spatial recaptures (Schmidt et al., 2022) and assessing for extreme distances moved across detected individuals to prevent unusual movements from biasing the σ estimate (Kendall et al., 2019).

We fit models in the software environment R (R Core Team, 2023), conducting model selection using the package oSCR (Sutherland et al., 2018). For each parameter, we fit univariate models to ensure variables captured the intended ecological effect, compare similar covariates (e.g., two hypotheses describing temperature), and understand univariate effect size and direction. We calculated pairwise correlations and retained only variables with correlations < 0.7 (Dormann et al., 2013). Because density was our primary interest, we considered detection as nuisance parameters and first identified best models for p ₀ and σ, using those in all hypothesized density models. For each model component (p ₀, σ, and density), we compared all possible combinations of the reduced variable set using Akaike’s Information Criteria (AIC, Burnham and Anderson, 2002).

We assessed constant, sex-specific, year-specific, and the combination of sex and year-specific variation in σ. We evaluated several hypotheses for influences on baseline detection: trap type and observer, plus trap covariates reflecting temperature in days before sampling and habitat conditions. For trap type, we competed hypotheses for categorizing trap type based on the initial reason for the stop (bighorn sign seen, targeted search location, bighorn seen, capture location, and opportunistic) comparing 2–5 categories. We tested the effect of the presence of two lead observers. We assessed temperature across the previous three and five days because warm days could drive bighorn to move closer to the river and increase detection. Habitat covariates summarized the area within 500 m of the trap site that we tested and included mean slope, maximum slope, terrain ruggedness, measured as the standard deviation of curvature (Ironside et al., 2018), exponentiated terrain ruggedness to account for neighboring steep cliffs, and summaries of habitat selection around the trap.

Because we expected more bighorn activity centers and higher detection of bighorn in places that they used more frequently, we developed a resource selection function (RSF) based on GPS locations from September, the month of pellet sampling, to use for detection and density covariates in the SCR model. Because this coincides with the breeding season when bighorn occur in mixed-sex groups and generally use the same areas, we included both sexes in a single model. We evaluated selection within a boundary that included the park plus the length of the Colorado River from Lees Ferry to the Diamond Creek takeout. Within that study area, we assigned each 30 m pixel as either used or unused depending on the presence of a telemetry point. We randomly assigned each of the unused sites proportionally to an individual. Thus, every location in the study area was included to completely describe the habitat in the analysis. We fit logistic models with a random intercept for individual, weighting each location based on the probability of fix acquisition, a spatially explicit probability based on topography that mitigates for the low fix success of some individuals (Ironside et al., 2018; Graves et al., 2024). We calculated aspect [northness = cos(aspect), eastness = sin(aspect)], the standard deviation of curvature and its quadratic, elevation and its quadratic, and solar radiation during the month of September as well its quadratic. We defined additional covariates based on the distance to a feature, namely distances to escape terrain >40° slope, to helicopter overflight paths, to trails, and to the river. For each distance covariate we also assessed a model with the effect of distance log-transformed, which can effectively model a threshold effect. After removing correlated variables, we fit univariate models. We eliminated the distance to helicopter covariate because the coefficient suggested attraction to helicopters, which is unlikely (Bleich et al., 1994), suggesting this was a spurious correlation (Wisdom et al., 2020). We fit a global model with all remaining variables. To use the RSF covariate as covariates on detection in the SCR model, we extracted mean and maximum values of predicted selection for each trap (location where boat stopped), competing linear and quadratic forms of RSF variables.

For the density portion of the model, we evaluated whether the mean or max of September RSF in a 500m radius, distance to the river, or both would best constrain activity centers to areas used by bighorn. We used binary variables to assess whether density changed inside the park boundary or in areas where domestic sheep grazing was possible based on land ownership. We considered three temporal model forms for changes in density: no time dependent covariate, linear yearly change, and year-specific density. We also assessed sex-specific variation in density.

Following model selection, we refit the best model from the maximum likelihood approach using a Bayesian approach to get credible intervals for estimated abundance (Woodruff et al., 2021), as detailed in Supplementary Appendix 1 . Abundance is a derived parameter based on spatially explicit density estimates. This model did not have separate intercepts for density by year, nor did it incorporate a trend term for density, as those model structures were less supported in model selection (see Results). The Bayesian model combined 1) delineation of individual activity centers across an inhomogeneous state space with covariates modeling the spatial variation in population density, 2) an observation model with categorical covariates to allow for variation in detection probabilities across J traps and K sessions, and 3) data augmentation to model the true population size. To implement data augmentation, we included an excess of all-zero encounter histories for individuals that were not detected during sampling but could be present in the population, M=1,200 individuals (Royle and Young, 2008; Royle et al., 2018; Woodruff et al., 2021). Our formulation assumed activity centers, s, for each individual, i, did not change over these 5 years. It also assumed the relationship of covariates with density was constant over time. We fit the Bayesian model using NIMBLE in the statistical computing environment R v.4.3.0 (de Valpine et al., 2017, de Valpine et al., 2023; R Core Team, 2023) and carried out 5 chains of MCMC sampling over 20,000 iterations with no thinning, of which we excluded the first 5,000 iterations for adaptation and burn-in.

From 2010–2018, we summarized 1,047 unique observations of bighorn groups, totaling 3,176 sightings of individuals, recorded on 94 river trips (private, commercial, and bighorn survey) to describe bighorn distribution along the river ( Figure 2 ). Bighorn counts averaged ~30 individuals per trip and were recorded between river miles 1 and 261. Overall patterns of distribution during this period appeared the same as the longer time series including pre-study dates (2000–2018). Our dedicated seven bighorn survey trips along the river (RM 0–RM 226) recorded 165 unique observations of bighorn groups totaling 535 bighorn sightings, and under the assumption that bighorn did not move down-canyon faster than our boats, yielded minimum counts that declined from 145 individuals (2011) to a low of 56–61 individuals (2012–2014), before increasing to 88–92 bighorn in 2015–2016 ( Table 1 ). The number of individuals detected via direct genetic identification of individuals from fecal pellets declined more slowly, from 134 in 2011 to 65 in 2013, before a rapid increase in genetic detections to 128 individuals in 2014 and nearly 200 in 2015, although the trap type varied substantially across years as did the numbers of individuals detected by trap type. For example, the number of traps where the reason for stopping involved seeing either bighorn or other sign varied from 17 to 49 across years. More stops in 2015 resulted in more samples collected. In contrast, fall lamb–ewe ratios generally declined from 50:100 in 2012 to an average of 24.3:100 from 2014–2016 ( Table 1 ).

We observed bighorn with signs of disease beginning in 2012, when a collared ewe died at RM 6 on RR and a necropsy revealed bronchopneumonia, although a ewe captured in 2010 exhibited nasal discharge. Reports of animals with symptoms congruent with respiratory disease (hereafter “sick”; see Methods) increased after 2012 and the combination of observations of symptoms and diagnostic samples suggested that a widespread disease event occurred with the maximum number of sick (n = 11) and dead (n = 16) bighorn reports in 2014, ranging from RM 8–196 on RR and RM 79–216 on RL ( Figure 3 ). The number of coughing and dead bighorn increased each year between 2012 and 2014 ( Supplementary Tables A1, A2 ). From 2011–2016, non-native respiratory disease pathogens were detected in 77% of dead and captured bighorn sampled for disease diagnostics ( Table 2 ). Testing of samples from collared bighorn indicated that 81% had been exposed to M. ovipneumoniae (serology) and 56% had active M. ovipneumoniae infections (PCR; Table 2 ). Most bighorn mortalities reported were not possible to investigate due to remoteness. Seven post-mortem investigations were completed, of which 2 yielded viable samples that both tested positive for M. ovipneumoniae, but most carcasses examined were too heavily scavenged or autolyzed to determine cause of death or test for M. ovipneumoniae or other pathogens ( Supplementary Table A2 ). Strain typing for M. ovipneumoniae, conducted by WADDL on 3 samples from collared animals at 4 loci (16S, gyrB, rpoB, IGS), indicated that 2 strains were detected in Grand Canyon during this study, none of which has been detected outside Grand Canyon (WADDL, unpublished data).

Between October 2010 and September 2016, we captured 25 individual adult bighorn (16 female and 9 males) between RM 6 and 214. Mean age across all sexes was 5.6 years. Estimated age of females at capture ranged from 3.5–7.5 years old (n = 16, mean = 5.5 years old) and males ranged from 2.5–8.5 years old (n = 9, mean = 5.9 years old). We captured 14 ewes and 5 rams on RR, and 2 ewes and 4 rams on RL ( Figure 4 ). All captured bighorn were collared except one ewe on RR. Two collars did not acquire GPS fixes and were not included in the RSF analysis but were monitored via VHF for 627–832 days. Collars with functioning GPS (n = 22) averaged 480 days (range = 178–754 days) and accumulated 10,560 collar-days (one collar day = one bighorn wearing a collar for one day) with 25,700 GPS fixes ( Supplementary Figure A1 ; Graves et al., 2024). Overall GPS fix success rates were poor, averaging 0.33 of programmed acquisitions (range = 0.07–0.72). Bighorn typically inhabit regions in the park with complex rugged topography and narrow sky views that constrain GPS fix success rates (Ironside et al., 2017). Although GPS fixes on individual collared bighorn were occasionally acquired on the opposite side of the river, these were single location events, and determined to be associated with GPS error based on topography, poor sky view, and pre- and post-locations.

For the global September RSF model, distance to the river, distance to escape terrain, and northness had the largest effects, with all covariates statistically significant ( Table 3 ). Bighorn had higher use of areas closer to the river closer to escape terrain, on southerly and easterly aspects, and in areas with moderately low variation in curvature and medium amounts of solar radiation ( Supplementary Figures A2–A4 ). Nearly 60% of September GPS locations were <500 m from the river, ~83% were <1 km from the river, and >90% were <1.5 km from the river.

For the population estimation dataset, we collected 1,539 fecal pellet samples ( Supplementary Figure A5 ; Graves et al., 2024). We successfully genotyped 1,250 samples at sufficient loci (4 minimum; Table 4 ) to resolve individual identity, yielding 453 individuals ( Table 1 ). After merging that dataset with the Creech et al. (2017) dataset (n = 262 individual genotypes, of which 208 were already represented within the population estimation dataset), the final combined dataset included 529 individuals, of which 257 individuals were genotyped at up to 16 loci (yielding 224 individuals with at least 15 loci resolved) and 273 were genotyped at up to 9 loci (yielding 233 individuals with 9 loci resolved); the full distribution of locus number and genotypes used for estimation of genetic structure is presented in Supplementary Figure A6 .

Based on individual identification from genotypes, we detected only 6 (3 male, 3 female) bighorn on both sides of the Colorado River, all between RM 136–141 ( Figure 4 ). All other redetections occurred on the same sides of the river and distances between individual redetections across the population estimation dataset ranged from 0–38 km ( Table 5 ; Supplementary Figure A7 ).

Estimates of the number of genetic clusters (K) using STRUCTURE with no informed priors on population or sampling location suggested K = 6 best represented structure within the data. At K = 6, the maximum q values ranged 0.97–0.98; at K = 7, the maximum q values per cluster ranged 0.77–0.98. Thus, applying the high assignment criterion, K = 6 was the highest supported value of K. ΔK values ( Supplementary Figure A8 ) suggested divisions within the data primarily at K = 2, 3, and 6. At K = 2, the primary division in the data reflected river side ( Supplementary Figure A9 ). After spatially defining populations based on individual assignments at K = 6 ( Figure 5 ), tests for Hardy–Weinberg proportions indicated that with Bonferroni correction for 16 tests, only one locus from the 16-locus data showed consistent violations. That locus, BL4, showed evidence of selection in another bighorn study area (Epps et al., 2018), but we retained it in this analysis. Of the six clusters, described hereafter as subpopulations, only 1 subpopulation (RL-East) showed 2 loci departing from HW proportions after Bonferroni correction for 16 tests; all other subpopulations showed 0 to 1 locus departing from HW proportions, suggesting that subpopulations we defined did not contain substantial additional substructure. Population genetic distance (F _ST) ranged from 0.022, between RR-Mid and RR-West subpopulations, to 0.126, between RL-East and RR-West ( Table 6 ).

Tests using STRUCTURE indicated migrants from both current and previous generations were detected among subpopulations on the same side of the river as well as among those separated by the river, particularly in the central section ( Supplementary Table A3 ). Estimation of migration rates among subpopulations in the previous generation using BIMr best supported a model where both factors (distance along the river and the river itself) limited migration. The run with the lowest Bayesian deviance estimated the posterior probability of the two-factor model at 0.75, compared to 0.10 for the null model with neither factor, 0.05 for the river alone, and 0.11 for distance alone; the model with river and distance was strongly supported in all 10 replicate simulations (mean posterior probability of the river and distance model was 0.72, range: 0.43–0.92). Migration rates (estimated for previous generation) were highest among RR-Mid and RR-West subpopulations on the same side of the river, but a high migration rate also was observed across the river between the middle subpopulations ( Table 6 ). Self-migration estimates clearly showed that both RR-Mid and RL-Mid subpopulations experienced a relatively high rate of migration (in the genetic sense) to and from other subpopulations ( Table 6 ).

Spatial recaptures ranged from 14 to 35 individuals per year ( Supplementary Table A4 ), suggesting sampling was sufficient to avoid bias (Schmidt et al., 2022), as sample size can influence bias and precision of population estimates from SCR approaches (Efford and Boulanger, 2019; Dupont et al., 2021). Because we used a single sampling session in each year, recaptures within sites did not exist in the analysis. We excluded one detection of an individual bighorn that was 38 km away from another detection ( Supplementary Figure A7 ; Table 5 ), exceeding 5 standard deviations of the maximum distance, because rare extreme movements can bias estimates of σ (Kendall et al., 2019). We retained movements <20km in the same year because we detected the same animals at this distance across different years for multiple animals, suggesting that this was a regular movement pattern. The proportion of males detected genetically in the SCR dataset ranged from 39–61% ( Table 1 ). Recaptures across years, albeit not used analytically, were generally consistent, with 23–39% of individuals detected 2012–2015 previously detected ( Supplementary Table A5 ).

Two models had high parsimony (<2 ΔAIC, next best model ΔAIC = 4.5), sharing identical form except that density either was constant (best model) or declined linearly with year. In the Bayesian fit of the best model, estimates generally aligned with likelihood-based estimates. Mean annual point estimates of abundance ranged from 536–552 with 95% credible intervals ranging from 451–647 across years ( Table 1 ). Both top models included the maximum September RSF value, distance to river, and park boundary indicator. More animals occurred in areas of better habitat, inside the park, and near the river, consistent with bighorn observations ( Table 7 ; Figure 6 ). Baseline detection included trap type as 4 categories and presence of 2 experienced observers. Detection was highest at sites where the team stopped because they saw bighorn or bighorn sign, followed by predefined search areas (2014–2015), systematic sites searched based on a congruence of good habitat and the ability to stop the raft (2011–2013), and finally by opportunistic locations. The best model for σ was session-specific, with the lowest estimated activity range in 2013 ( Table 7 ). Sex differences in density or detection were not supported.