Sex‐specific population dynamics of ocelots in Belize using open population spatial capture–recapture

Abstract

We used open population, spatial capture–recapture (SCR) models to estimate sex‐specific density, survival, per capita recruitment, and population growth rate of ocelots (Leopardus pardalis) at five sites in Belize with up to 12 yr of data per site. Open population SCR models enabled us to separate survival and recruitment from migration using an ecologically realistic, spatially explicit movement model. Yearly survival probability across 4 broadleaf forest sites was estimated at 0.73–0.84 for males and 0.81–0.87 for females, with no clear indication of sex differences. Yearly per capita recruitment was estimated across four broadleaf forest sites at 0.06–0.08 recruits/N for males and 0.09–0.12 recruits/N for females, again with no clear indication of sex differences. At a pine forest site with a population comprised largely of males, survival and recruitment estimates were similar to the broadleaf sites. Population densities in the broadleaf forest sites ranged from 6.5 to 14.7 ocelots/100 km², and 0.9–2.5 ocelots/100 km² in the pine forest site, with strong evidence of a female‐biased sex ratio in the broadleaf sites and a male‐biased sex ratio in the pine forest site. We also found strong evidence that female within‐year space use at the broadleaf sites was smaller than that of males, and that within‐year space use at the pine forest site was larger than that at broadleaf sites. Between‐year home‐range relocation at broadleaf sites was of a similar spatial scale as within‐year space use, consistent with philopatry. We found evidence of a small population decline (posterior probability > 0.9) at two of four broadleaf sites; however, given the level of uncertainty about decline magnitudes, we suggest continued monitoring of these sites to increase site‐years and gain further precision on population growth rate estimates. Estimating demographic parameters at large spatial and temporal scales is important for establishing reliable baseline estimates for future comparison and for understanding changes in population dynamics. Long‐term data sets like those we collected are of particular importance for long‐lived species living at low densities and large spatial scales, where not many individuals are exposed to capture in any one year.

Introduction

The ocelot (Leopardus pardalis) is distributed across South and Central America, the islands of Trinidad and Margarita, Mexico, and southern Texas in the United States. The species is currently listed under the IUCN Red List of Threatened species as least concern and Appendix I of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES), which prohibits all international trade of skins and live animals. Retaliatory killing due to depredation on poultry (Amador‐Alcalá et al. ), along with habitat loss and fragmentation, is currently the main threat to ocelot survival (Sunquist and Sunquist ). In Belize, and other Neotropical countries, deforestation stemming from large‐scale agriculture and infrastructure development (Young , Nogueira and Nogueira‐Filho , Aide et al. ), as well as increased human populations, threaten the long‐term survival of ocelots. Knowledge of contemporary population trends across multiple locations is necessary to determine when and where conservation action may be needed. Further, an understanding of how survival and recruitment contribute to population growth rate is necessary to provide a baseline for comparison to future population trends. Of particular importance is the recruitment rate since ocelots naturally have low fecundity (e.g., low reproductive rates, long inter‐birth intervals, and small litter sizes).

Ocelots have been studied in multiple locations across their range; however, most studies focusing on population ecology are limited to estimation of abundance and density (Dillon and Kelly , , Martínez‐Hernández et al. , da Rocha et al. , Gómez‐Ramírez et al. , Satter et al. ). To date, however, only three studies have estimated parameters associated with ocelot population change through time (e.g., survival, recruitment, and population growth rate). For example, in Texas, survival was estimated in a single year via VHF monitoring (Haines et al. ) and in Mexico (over three years; Gómez‐Ramírez et al. ) and Belize (over six and seven years; Satter ) using camera‐trap capture–recapture data using robust design models (Kendall et al. ). However, no study has yet estimated an ocelot population growth rate or investigated how survival and per capita recruitment rates contribute to population trajectory through time, which we attribute to lack of long‐term, capture–recapture data sets that allow for enough statistical power to estimate these parameters with adequate precision.

While it is possible to estimate survival, per capita recruitment, and population growth rate with as few as two years of capture–recapture data, the precision of these parameter estimates for low‐density populations of K‐selected species will frequently not be of adequate precision to reliably guide conservation action (White ). Further, temporal autocorrelation in these population parameters may lead to parameter estimates over short time periods that are inaccurate estimates of long‐term averages due to a loss of statistical efficiency (Bence ). Longer‐term, capture–recapture surveys maximize the likelihood that population demographic estimates will be precise enough to reliably assess the current population trajectory and provide more accurate baseline estimates of survival and recruitment with sufficient precision to detect future deviations. Further, because population parameters may vary across space, estimating such parameters across multiple sites will provide a more thorough understanding of the degree of variability across sites within a region.

A second way to improve our understanding of a species’ population dynamics is to use more comprehensive statistical models. A perennial challenge in classical open population capture–recapture is the separation of the survival and emigration processes. Further, the estimation of recruitment in the presence of immigration can only be achieved using robust design models when there is no emigration (Cooch and White ) or when the population can be divided into age classes where recruitment from the juvenile to adult state can be observed (Nichols and Pollock ). Open population, spatial capture–recapture (SCR) models offer a spatially explicit, and thus more realistic, representation of the population dynamics processes that allow both survival and recruitment parameters to be estimated without information on age class (Gardner et al. , Chandler and Clark ; with caveats, see Discussion). This is achieved by modeling activity center relocation between primary periods, which is only possible using models that capitalize on the information contained in the spatial locations of capture events. A by‐product of this integrated movement model is that it allows for inference about the process of activity center relocation itself, and how they may differ across habitat types or between population subgroups.

Our objective was to use a unique long‐term, multi‐site data set to estimate ocelot density, survival, per capita recruitment, inter‐ and intra‐year activity center movement scales, and population growth rates over time (up to 12 yr) and space (five sites) in Belize with sex‐specific, open population SCR models. Four sites, La Milpa, Hill Bank, Cockscomb, and Gallon Jug, were comprised of broadleaf forest—considered high‐quality ocelot habitat—while Mountain Pine Ridge was comprised of higher elevation pine forest, considered marginal ocelot habitat supporting a lower population density (Satter et al. ). This data set was previously used to estimate site‐by‐year population density (Satter , Satter et al. ) using independent, closed SCR models, with a subset of the data previously used to estimate survival using robust design (Satter ). These previous site‐by‐year density estimates (Satter et al. ) suggested population growth rate, λ, was roughly stable in Gallon Jug and Hill Bank, but small increases or decreases in population size could not be ruled out. The year‐by‐site density estimates at La Milpa and Cockscomb were more ambiguous, with a downward trend, but high sampling variance, precluding determination of general population trajectory. Thus, one of our main objectives was to share information about density among years through a population growth model to determine the population trajectories at these five sites with greater precision.

Because ocelots are a K‐selected species, we expected high survival probabilities and correspondingly low per capita recruitment rates. In the presence of abundant, stable, and evenly distributed prey sources, females often have smaller exclusive territories than males (Sunquist and Sunquist ). Given that males compete with one another for space use and access to females, we expected any sex‐specific survival differences to be in the direction of higher survival for females. Then, if sex ratios were stable through time at each site, we expected a lower per capita recruitment rate among females to fill the fewer female vacancies due to higher adult survival. Finally, we expected male, within‐year space use in the broadleaf sites to be larger than females as previously documented (Satter et al. ) and between‐year activity center movement of a similar spatial scale as within‐year movement, corresponding to small between‐year activity center shifts and reflecting a general trend of philopatry.

The expected population dynamics patterns at Mountain Pine Ridge were less clear because it is largely comprised of the seemingly non‐preferred pine forest habitat as evidenced by the very low number of ocelot detections, surrounded by the preferred broadleaf forest. But this habitat matrix may act as a temporary refugia for subadult males or other subordinate individuals who cannot obtain a territory exclusively in broadleaf forest. Therefore, we predicted that the sex ratio would be skewed toward males, and individuals would have larger within‐ and (especially) between‐year spatial scale parameters, and lower survival than at the broadleaf sites. Finally, we predicted the majority of individuals detected at this pine forest site would live on the periphery of the site, closer to, or largely within, the neighboring broadleaf forest.

Methods

Study area

Our data set comes from long‐term camera‐trap surveys at 5 sites in Belize, Central America, conducted from 2004 to 2016 (Fig. ). The surveys were conducted in Cockscomb Basin Wildlife Sanctuary (2004–2008, 2011–2015), Mountain Pine Ridge Forest Reserve (2004–2015), Gallon Jug Estate (2013–2016), and two sites within the Rio Bravo Conservation and Management Area: Hill Bank (2010–2015) and La Milpa (2008, 2010–2015). Cockscomb encompasses 425 km² broadleaf tropical moist rainforest in south‐central Belize, with elevations of 50–1120 m (Silver et al. ). Mountain Pine Ridge located in central‐west Belize is an ~434 km² pine (Pinus sp.) dominated system with some smaller areas of shrub and broadleaf moist forest, and elevations of 120–1017 m (Wultsch et al. ). Gallon Jug, located in northwestern Belize, is ~538 km² (including some parts of Yalbac Ranch and Laguna Seca) and is composed primarily of lowland broadleaf moist evergreen seasonal forests (Miller ), with an elevation of ~40–160 m. Finally, both La Milpa and Hill Bank are primarily broadleaf forest sites with an elevation ranging from 40 to 160 m. However, La Milpa has more upland broadleaf forest, and Hill Bank is lowland broadleaf forest with large areas intermixed with freshwater swamp and pine savanna. For a more detailed summary of all five study sites, see Satter et al. ().

Locations for ocelot camera‐trap surveys for five sites in Belize, Central America, 2004–2016. Shapefiles were adapted from the Biodiversity and Environmental Resource Data System of Belize (BERDS).

Field methods

At each site, ~20–50 paired camera‐trap stations were established (depending on year and site), with sampling periods no longer than three months. Based on home‐range estimates for ocelots, pumas, and jaguars, camera stations were spaced systematically between 2 and 3 km² apart (Emmons , Silver et al. , Di Bitetti et al. , Haines et al. , b), on opposite sides of trails and logging/or multi‐use roads to increase the likelihood of photographing individuals on both flanks (Dillon and Kelly ). We identified individual ocelots by their distinct coat patterns, differentiated sex by a presence or absence of testicles, and recorded sex as unknown if it could not be determined for any particular individual (see Satter et al. for further details). To increase sample size when constructing the capture histories, we appended captures of individuals photographed on only one flank (the larger of the left‐only or right‐only data set) to captures of individuals simultaneously photographed on both flanks at least once. While this practice has been shown to introduce low‐to‐moderate individual heterogeneity in capture probability (Augustine et al. ), the magnitude of the effect is small in our data set because multiple years substantially increased our probability of obtaining at least one simultaneous capture of both flanks for each individual; thus, we had relatively few single‐sided only photographs.

Statistical methods

We used an open population, spatial capture–recapture model in a Bayesian framework (Gardner et al. , Chandler and Clark ) to estimate combined and sex‐specific ocelot population parameters (survival, per capita recruitment, and population growth rate) separately for each of the five sites. At each site, the field methods produced robust design structured data sets where each year was a primary period and each day (of the 2‐ to 3‐month survey period) within year was a secondary period. The secondary periods allowed for yearly closed population abundance to be estimated in the presence of imperfect detection, while the primary periods allowed for population parameters to be estimated in the presence of uncertainty in abundance and possible immigration and emigration. To separate emigration and immigration from survival and recruitment, we allowed for population redistribution between primary periods (Gardner et al. , Ergon and Gardner ). At each site, the camera‐trap detection history was formatted in a three‐dimensional array of dimension n × J^max × L, where n is the total number of individuals captured over all years at a site, J^max is the maximum number of traps in any year at a site, and L is the number of years the site was sampled. We used the R package OpenPopSCR (Augustine ) to estimate the parameters of the open population SCR model from these data sets, which we describe in more detail below.

Process model

We used a density‐independent population growth model, modified from Chandler and Clark () to allow parameters to vary by sex. The population parameters of the model are φ^m and φ^f, the male and female yearly survival probabilities, and γ^m and γ^f, the male and female yearly per capita recruitment rate. These yearly, sex‐specific per capita recruitment rate parameters are the number of individuals of each sex added per year per total abundance, for example, γ^m is the number of males added per N each year. We denote total, male, and female abundance in year l as N_l, $N_l^m$, and $N_l^f$, respectively. Then, starting in year two, the male and female abundance in each year is determined by their abundance in the previous year following $N_{l+1}^m={\mathrm{\phi }}^m{N}_l^m+{\mathrm{\gamma }}^m{N}_l$ and $N_{l+1}^f={\mathrm{\phi }}^f{N}_l^f+{\mathrm{\gamma }}^f{N}_l,$ and the total abundance in each year is the sum of the sex‐specific abundances. Associated with each individual i in year l is a two‐dimensional activity center, s_il, with the year specificity allowing each individual to relocate between years. The matrix S is thus an array of dimension N × L × 2 storing the X and Y coordinates of each individual in each year. We make the typical assumption of spatial uniformity of activity centers (Royle et al. ) in the first year, that is, $s_{i,1}\sim \text{Uniform}\left(\mathcal{S}\right),$ where $\mathcal{S}$ is a continuous two‐dimensional state space representing the area in which the population lives, which we defined by buffering the maximal extent of the trapping arrays in the X and Y dimension at each site pooled across years by at least three times the estimated detection function spatial scale parameter, σ (Royle et al. ). This rectangular state‐space definition created state spaces of areas 2010.9, 1526.2, 1333.4, 925.0, 1100.3 km² for Cockscomb, Mountain Pine Ridge, Gallon Jug, La Milpa, and Hill Bank, respectively.

Conditional on the initial activity center locations in the first year, the activity center of individual i in year l + 1 depends on its location in year l following $s_{i,l+1}\sim \text{Bivariate Normal}\left(s_{i,l},{\mathrm{\sigma }}^{L}I\right)$ truncated by the state‐space extent (Royle et al. ). Note, the year‐level spatial scale parameter, which determines how far individuals can relocate, is denoted σ^L to distinguish it from the detection function spatial scale parameters defined below. Due to poor mixing of this parameter with typically sparse capture–recapture data sets, we did not consider that it could vary by sex. Also associated with each individual is its sex, stored in the N_T‐dimensional vector C, where N_T is the total number of individuals alive in the population across all years. We assume the sex of individual i is distributed as c_i ∼ Bernoulli(p^sex),  where a value of 0 indicates a male, a value of 1 indicates a female, and p_sex is the probability any randomly selected individual in the superpopulation, N_T, is female. The expected sex ratio in each year can then be derived from the process model parameters, and the realized sex ratio in year l is a derived parameter, $p_l^{R.\mathrm{sex}}=N_l^f/N_l$. See Royle et al. () for a likelihood‐based treatment of class‐structured SCR models.

We also considered a simplified version of the process model described above, using a single parameter for both per capita recruitment and survival to remove the sex specificity. We considered this model for two reasons. First, this model was more appropriate for the Mountain Pine Ridge site with very sparse data consisting almost exclusively of males. Second, this model provided more precise parameter estimates for per capita recruitment and survival for the other four sites where strong evidence of sex specificity in these parameters was not observed.

Observation model

In each year l, ocelots were observed at camera locations X_l over K_l occasions, recorded as binary detection events, and then summed over occasions. Therefore, we assume the individual by trap‐by‐year number of observations was distributed as $y_{\mathit{ijl}}^m\sim \text{Binomial}\left(p_{\mathit{ijl}}^m,K_l\right)$ and $y_{\mathit{ijl}}^f\sim \text{Binomial}\left(p_{\mathit{ijl}}^f,K_l\right)$ for males and females, respectively, where $p_{\mathit{ijl}}^m$ and $p_{\mathit{ijl}}^f$ are the sex‐specific individual by trap‐by‐year capture probabilities, which we assume are determined by a hazard half normal detection function conditional on the yearly activity centers (Royle et al. ). This detection function has parameters λ₀, the baseline detection rate, and σ, the spatial scale parameter. The sex‐specific expected number of counts for individual i at trap location j in year l on each occasion was then ${\mathrm{\lambda }}_{\mathit{ijl}}^m={\mathrm{\lambda }}_0^{m}exp\left({\left|\left|{\mathbf{s}}_{\mathit{il}}-{\mathbf{x}}_{\mathit{jl}}\right|\right|}^2/2{\left({\mathrm{\sigma }}^m\right)}^2\right)$ and ${\mathrm{\lambda }}_{\mathit{ijl}}^f={\mathrm{\lambda }}_0^{f}exp\left({\left|\left|{\mathbf{s}}_{\mathit{il}}-{\mathbf{x}}_{\mathit{jl}}\right|\right|}^2/2{\left({\mathrm{\sigma }}^f\right)}^2\right)$ for males and females, respectively. Then, the sex‐specific expected number of counts were transformed to sex‐specific capture probabilities following p_ijl = 1 – exp(−λ_ijl). The only exception to the observation model described above is that the detection function parameters at Mountain Pine Ridge were not assumed to be sex‐specific due to sparse data and a predominantly male data set. The observed sexes were stored in vector C, now of length n_L, the total number of individuals captured across all years, with the possibility that the sex of some individuals is not observed. Finally, if a specific site was not surveyed in a specific year (2009–2010 at Cockscomb and 2009 at La Milpa), no data from that year contributed to the estimation of parameters, and the abundances in those years are estimated through the temporal dependence in the process model (see Chandler and Clark for previous application). To avoid confusion, we note here that the symbol λ is used for both the baseline detection rate and population growth rate to maintain consistency with the literature, with the former accompanied by a zero subscript: λ₀.

Inference via MCMC

We estimated the model parameters with typical SCR, MCMC algorithms using data augmentation (Royle et al. ) via the OpenPopSCR R package (Augustine ). At each site, the three‐dimensional capture history was augmented up to dimension M × J^max × L, where M is chosen to be much larger than the expected number of individuals alive across all years at a site. The observed sex vector, C, was also augmented up to length M, which is used to estimate the sexes of uncaptured individuals and captured individuals of unknown sex. In open population SCR, an M × L indicator matrix z is used to indicate whether individual i is in the population in year l (Gardner et al. , Chandler and Clark ). The yearly realized abundances are parameters derived from z by summing its columns, producing samples from the posteriors of yearly realized abundance. We used a sex‐specific version of this algorithm where the posterior samples of yearly sex‐specific realized abundance were derived following $N_l^m={\sum }_i{z}_{\mathit{il}}\times I\left(c_i=0\right)$ and $N_l^f={\sum }_i{z}_{\mathit{il}}\times I\left(c_i=0\right)$, where I() is an indicator function evaluating whether the sex of individual i is male (0) or female (1). Combined and sex‐specific yearly densities were then derived by dividing the appropriate abundance by the area of the state space. We refer readers to Chandler and Clark () for how survival and per capita recruitment are estimated using the z matrix. Combined and sex‐specific realized yearly population growth rates were derived from the sex‐specific abundances following λ_l = N_l/N_l−1, ${\mathrm{\lambda }}_l^m=N_l^m/N_{l-1}^m$, and ${\mathrm{\lambda }}_l^f=N_l^f/N_{l-1}^f$. To pool these yearly realized population growth rates and realized sex ratios into single, more precise estimates for each site, we used an inverse variance weighted random effects estimator (Borenstein et al. ) commonly used for meta‐analyses after determining that there was no evidence of temporal correlation in these parameters across years. All derived parameters (realized density, realized population growth rate, and realized sex ratio) came from the model with sex‐specific population parameters, except for the Mountain Pine Ridge site. All point estimates were obtained using the posterior mode, and all 95% credible intervals were obtained using the 95% highest posterior density (HPD) intervals.

Mountain pine ridge realized density plot

In order to test the hypothesis that individuals exposed to capture at the pine forest site predominately lived in or near the surrounding broadleaf forest, we produced spatially explicit, posterior density plots, split into four to three‐year periods to better visualize the change in realized density in and near the pine forest through time. Royle et al. () demonstrate how the activity centers and data augmentation vector z can be used to estimate the abundance within any region by summing the subset of the z_i vector that is located within the prescribed region. Because we wanted to depict density variation across the state space in a continuous manner, we instead produced kernel density plots of the activity centers for which z_il = 1 during each of the four to three‐year periods using a bandwidth of 2.5 km.

Results

Field results

Across all five sites and years, we had a total sample size of 74,854 trap nights and 2203 ocelot detections resulting in a trap success rate of 2.94 detections per 100 trap nights over all sites combined. Trap nights at each site in any given year ranged from 815 to 2894 and the cumulative number of ocelot detections at each site across years ranged from 111 to 813. The number of unique ocelots captured at each site in any given sampling year varied from 1 to 51, with spatial recaptures (i.e., captured at >1 camera station) ranging from 0 to 24 (Appendix S1: Table S1). In addition, across all sites and years, we identified a total of 322 adult ocelots: 148 males (M), 148 females (F), and 26 of unknown sex (UK) across all years in Cockscomb (39M:35F:11UK), Hill Bank (19M:21F:1UK), La Milpa (41M:38F:7UK), Gallon Jug (31M:50F:3UK), and Mountain Pine Ridge (18M:4F:4UK).

Statistical results

Yearly per capita recruitment point estimates ranged from 0.064 to 0.080 for males and 0.085 to 0.116 for females across the 4 broadleaf forest sites (Fig. ; Appendix S2: Table S1). The estimates without sex specificity, including the pine forest site ranged from 0.070 to 0.101. Judging by the overlap of 95% HPD intervals, there was not a strong indication that per capita recruitment varied between the sexes or across sites, although all 4 female point estimates were larger than those for the males (Fig. ). Yearly survival probability point estimates ranged from 0.732 to 0.837 for males and 0.809 to 0.868 for females across the four broadleaf forest sites (Fig. ; Appendix S2: Table S1). The estimates without sex specificity including the pine forest site ranged from 0.767 to 0.832. There was also no strong indication that yearly survival varied between the sexes or across sites. The precision of these per capita recruitment and survival parameter estimates as measured by the 95% HPD width varied as a function of both the number of years considered (surveyed years and un‐surveyed years between the surveyed years) and the population density. Estimates were most precise at Cockscomb with a period of 12 yr and then La Milpa with a period of 8 yr. Gallon Jug estimates were the next most precise, despite a shorter sampling duration than Hill Bank (four vs. six years), due to the substantially higher density, and thus, less sparse yearly data sets, at Gallon Jug. Mountain Pine Ridge had the least precise estimates due to the very low population density and the resulting sparse yearly data sets, despite having a length of 12 yr.

Yearly combined and sex‐specific per capita recruitment estimates (posterior modes and 95% highest posterior density intervals) at each site in Belize. These parameters are the number of males and females added each year per total N of both sexes. In the sex‐specific model, these parameters vary by sex, and in the combined model, the same number of males and females is recruited into the population each year. Abbreviations are CC, Cockscomb; HB, Hill Bank; LM, La Milpa; GJ, Gallon Jug; MPR, Mountain Pine Ridge.

Yearly combined and sex‐specific survival estimates (posterior modes and 95% highest posterior density intervals) at each site in Belize. These parameters are the probability that an individual survives 1 yr. Abbreviations are CC, Cockscomb; HB, Hill Bank; LM, La Milpa; GJ, Gallon Jug; MPR, Mountain Pine Ridge.

The inverse variance weighted realized population growth rate point estimates (Fig. ; Appendix S2: Table S2) ranged from 0.89 to 0.99 for males and from 0.96 to 1.09 for females across sites. The point estimates for the combined sex population growth rate ranged from 0.91 to 1.05. There was not a strong indication that population growth rates varied between the sexes or across sites. The posterior probability that the total population size was declining (λ < 1) was 0.90, 0.44, 0.95, 0.19, and 0.99 for Cockscomb, Hill Bank, La Milpa, Gallon Jug, and Mountain Pine Ridge, respectively. The posterior probabilities that these populations were declining by more than 5%/yr (λ < 0.95) were 0.14, 0.11, 0.40, 0.06, 0.85, respectively. The precision of the inverse variance weighted population growth rate estimates as measured by the 95% credible interval width varied as a function of the number of years considered, except for Mountain Pine Ridge, where density was very low and data were sparse (Fig. ). The yearly realized λ estimates used to compute the inverse variance weighted estimates can be seen in Appendix S2: Tables S6 and S7 and Fig. S1.

Combined and sex‐specific realized population growth rate estimates (posterior modes and 95% highest posterior density intervals) for each site in Belize, obtained using an inverse variance random effects estimator on the yearly estimates at each site. λ = 1 indicates a stable population size through time. CC, Cockscomb; HB, Hill Bank; LM, La Milpa; GJ, Gallon Jug; MPR, Mountain Pine Ridge.

Combined realized population density point estimates ranged from 6.5 to 14.1 individuals per 100 km² across years at the broadleaf forest sites and from 0.9 to 2.5 in the pine forest site (Fig. ; Appendix S2: Tables S4 and S5). These point estimates declined through time at Cockscomb, Mountain Pine Ridge, La Milpa, and Hill Bank and increased through time at Gallon Jug; however, there was considerable overlap of the 95% HPD intervals between the first and last years surveyed at all sites and inference about population trends should be made from the inverse variance weighted population growth rate estimates. The female density point estimates tended to be higher than males at the broadleaf forest sites and lower than males at the high elevation pine forest site. The yearly realized sex ratio point estimates (Appendix S2: Table S8, Fig. S2) indicated a female bias in all years at the broadleaf sites and a male bias in all years at the pine forest site. The 95% HPD intervals on sex ratio did not overlap 0.5 (i.e., equal sex ratio) in any of the 4 yr at Gallon Jug, 4/6 yr at Hill Bank, and 8/12 yr at Mountain Pine Ridge, indicating female bias in the broadleaf forest sites and male bias in the pine forest site. The inverse variance realized sex ratio estimates (Appendix S2: Table S9) indicated female bias at all broadleaf sites (posterior probability P(female) > 0.5; Cockscomb 0.999, Hill Bank 0.999, Gallon Jug 0.997, La Milpa 0.998) and male bias at the pine forest site (posterior probability P(female) <0.5 = 1.000).

Yearly realized density estimates (posterior modes and 95% highest posterior density intervals) for each site in Belize, across all years monitored between 2004 and 2016. Note all Y‐axes are of the same height except for Mountain Pine Ridge.

The within‐year detection function spatial scale, σ, point estimates for the broadleaf sites ranged from 1.78 to 2.22 km for males and from 1.14 to 1.73 km for females across years (Fig. ; Appendix S2: Table S3). The combined sex σ point estimate for the high elevation pine forest site was 3.17 with a 95% credible interval that did not overlap any of the estimates from the broadleaf sites. The male and female σ 95% credible intervals did not overlap at any of the broadleaf sites except for Gallon Jug where the overlap was slight. The pooled‐sex activity center relocation spatial scale parameter, σ^L, estimates were similar for the broadleaf forest sites and the estimate for the pine forest site was larger than those of the broadleaf forest sites. The implied 1‐dimensional relocation kernels can be found in Fig. . The 5%, 50%, and 95% quantiles of movement distances between years based on the estimated bivariate normal redistribution kernels for Gallon Jug were 0.46, 1.74, and 3.64 km and 1.38, 5.08, and 10.57 km for Mountain Pine Ridge. The baseline detection rate estimates, λ₀, can be found in Table S3 of Appendix S2.

Within‐ and between‐year spatial scale parameters for ocelots in Belize. Parameters denoted DF are the within‐year detection function σ parameters, which were sex‐specific for all sites except for Mountain Pine Ridge. The between‐year dispersal parameters, σL (noted as Dispersal Combined), were not sex‐specific. CC, Cockscomb; HB, Hill Bank; LM, La Milpa; GJ, Gallon Jug; MPR, Mountain Pine Ridge.

The estimated dispersal kernel for ocelots implied by the σL estimates for Gallon Jug (representative of the lowland broadleaf forest sites), and Mountain Pine Ridge (the high elevation, low population density pine forest site) in Belize. Point estimates are indicated with solid lines, and 95% credible limits are indicated with dotted lines.

Finally, the posterior realized density at Mountain Pine Ridge broken into four time periods suggested that the activity centers of individuals living in this pine forest–broadleaf habitat matrix were predominantly in the broadleaf forest or at the edge of the pine forest near the broadleaf forest (Fig. ). Individuals with activity centers in the pine forest tended to live closer to streams. There was a decreasing posterior density within the pine forest from a high in the 2004–2006 period to a low in the 2010–2012 period, and a slight increase in the 2013–2015 period.

Spatially explicit, realized population density (summed activity center posteriors) for ocelots at Mountain Pine Ridge, Belize, in three‐year intervals showing a declining density in the pine forest through 2012, with an increase during the next three‐year period. Darker blue indicates greater population density, traps (black dots) are weighted by the number of captures during each three‐year period, pine forest is delineated with a thick solid line, and the area south of the thin dashed line is contiguous broadleaf forest, while the area north of this boundary is fragmented broadleaf forest.

Discussion

We applied open population SCR models to an ocelot camera‐trapping data set spanning up to 12 yr across five sites in Belize for a total of 42 site‐years, 39 of which were surveyed. This large spatial and temporal data set allowed us to estimate population demographic parameters with greater precision than shorter‐duration studies and the multiple sites allowed for a more robust picture of ocelot population status in Belize than can be ascertained from the typical single‐site studies. The use of open population SCR models allowed us to separate survival from emigration and recruitment from immigration (Ergon and Gardner , Schaub and Royle ) in an ecologically realistic, spatially explicit manner, and allowed us to compare the spatial scale of movement and the scale of activity center relocation at sites. Further, the sex specificity allowed us to compare detection and population parameters between sexes.

As expected for a K‐selected species, estimated survival probabilities were high and estimated recruitment rates were correspondingly low, indicating little population turnover each year. Per capita recruitment was estimated at ≤0.1 recruits per N across sites, which corresponds to fewer than 10 recruits of each sex (fewer than 20 total) per year in a hypothetical population of size 100. Estimated survival probabilities were 0.81–0.83 at the broadleaf sites and 0.76 at the pine forest site, corresponding to the death of about 20 individuals (out of 100) per year. Thus, there was an estimated turnover in about 1/5 of the population per year due to births and deaths. We did not detect any differences in per capita recruitment or survival rates between sexes or across sites; however, even with the statistical power provided by many years of data at each site, we likely still did not have enough power to detect ecologically meaningful differences in these parameters, especially in per capita recruitment. For example, the mean 95% HPD interval width for combined sex per capita recruitment was 0.086, which is wide enough to cover both 0.05 and 0.10, corresponding to an expected addition of five or 10 individuals of each sex per year or 10 or 20 total individuals in a population of size 100. This is an ecologically significant difference that could easily determine whether a population was increasing or decreasing through time. The survival probability was estimated with greater precision because it is not estimated per capita, and therefore, the uncertainty in N is not incorporated into the uncertainty in survival. Likewise, there is likely lower uncertainty in the estimated number of realized recruits, which can be calculated from the z matrix and not a function of N (we did not calculate or present these); however, per capita recruitment is generally of more interest and is required to make comparisons across populations of differing sizes and is useful for predicting population sizes into the future, for example, using population viability analyses.

To our knowledge, we provided the first per capita recruitment rate estimates for ocelots so we have no basis for comparison of these estimates; however, we can compare our survival estimates to those from previous studies. First, our survival estimates were substantially higher than those produced by Gómez‐Ramírez et al. (), who estimated apparent survival (0.64 vs. 0.77–0.83 in this study), assuming no emigration from the study area. This suggests that emigration accounts for a significant proportion of individuals disappearing from a site between years and/or survival was actually lower in this Mexican population. We suspect the confounding of emigration and survival explains the majority of the difference between our survival estimates and those from Gómez‐Ramírez et al. (). Second, our survival point estimates are slightly lower than those for resident ocelots in southern Texas (0.87), but higher than those of transients (0.57; Haines et al. ). Using camera‐trap capture–recapture, we could not differentiate between residents and transients, but the average survival estimates from Haines et al. (, b) weighted by the proportion of residents and transients in their population (0.78 vs. 0.22) is 0.80, which is very similar to our point estimates for ocelots without consideration of the residency status. Finally, Satter () estimated survival from the same data sets as we used for Hill Bank and La Milpa using classical robust design methods, allowing for temporary emigration and produced estimates that were both slightly lower than those we produced (0.83 vs. 0.79 for Hill Bank and 0.78 vs. 0.71 for La Milpa), although our interval estimates largely overlapped and were of similar precision. The higher point estimates in our model with spatially explicit dispersal dynamics are consistent with the simulation study of Horton and Letcher () that found the robust design model to underestimate true survival when simulating data from a population with spatially explicit dispersal; however, we cannot say with certainty that this mechanism explains the difference in our results.

We provided the first estimates of ocelot population growth rates; however, even with many years of data at each site, there remains a considerable level of uncertainty about the population trajectories. One disadvantage of using the Bayesian open population model utilizing data augmentation is that expected and realized population growth rates are year‐specific‐derived parameters and thus cannot be pooled across years or directly treated as random effects. The strategy we employed to pool information using an inverse variance weighted random effects meta‐analysis method (Borenstein et al. ) relies on the assumption of independence in the point estimates across years. Independence is more likely to be met using the realized, rather than expected population growth rates, because they may deviate from the deterministic population growth trajectory. In our case, the realized λ values varied around a stable mean and were effectively independent through time (Appendix S2: Fig. S2), but using the random effects estimator assuming independence is not likely to be a strategy that will be generally applicable to species that are growing or declining.

The inverse variance weighted realized population growth rates were more precise than the year‐by‐site estimates (Fig.  vs. Appendix S2: Fig. S1), especially for the sites with 8 or more years of data (Cockscomb, La Milpa, and Mountain Pine Ridge), but there was still substantial uncertainty about the trajectory of population size in the broadleaf forest sites. The posterior probability that the total population size was declining (λ < 1) was high in three populations—Mountain Pine Ridge (0.99), La Milpa (0.95), and Cockscomb (0.90); however, the only population with a 95% credible interval that did not overlap and a high posterior probability of a greater than 5% decline per year (λ < 0.95) was Mountain Pine Ridge (0.85).

The site‐by‐year realized density point estimates (Fig. ) agreed with the realized population growth point estimates, with density at Mountain Pine Ridge declining over 100% and moderate, but stabilizing declines at La Milpa and Cockscomb. Realized population density point estimates were stable through time at Hill Bank and stable to increasing at Gallon Jug. The evidence of population declines of <5% per year at the broadleaf sites of La Milpa and Cockscomb should be interpreted along with the still substantial estimated population densities at these sites of 9.8 and 6.8 individuals/100 km², respectively, in the final survey. Further, the substantial overlap in the interval estimates between the first and last years of surveys at each site reflects substantial uncertainty that the populations actually declined during the period surveyed. Finally, these population growth estimates may overstate evidence for decline if there was temporal variability in survival and recruitment, which cannot be modeled with precision with an average of 17 individuals exposed to capture at each site/year. This is a fundamental limitation of sampling K‐selected, low‐density populations with camera traps.

Although we did not find strong evidence of sex‐specific population dynamics, we did find evidence of sex bias in density and space use dynamics. The estimated sex ratios at all broadleaf forest sites were female‐biased in all years, and the estimated sex ratio at Mountain Pine Ridge was always male‐biased. Male within‐year space use was greater than females at the four broadleaf sites as measured by the estimated detection function spatial scale parameters. Our results, which demonstrated evidence of higher female densities due to female‐biased sex ratios, and lower within‐year space use by females, may be a consequence of typical felid social system dynamics. For example, females often have smaller home ranges than males and tend to be philopatric, while males disperse as sub‐adults, which then creates higher clusters of related females imbedded within home ranges (Sunquist and Sunquist ). In addition, higher female densities may be partially explained by common trends in mammalian populations where behavior (e.g., aggression), disease, genetics, and physiological factors tend to skew wildlife population sex ratios toward female (Mills ).

The spatially explicit, open population model allowed us to quantify the magnitude of between‐year activity center relocation and infer that on average, ocelots are largely philopatric in the broadleaf forest sites, but less so in the pine forest site. Comparing the detection function spatial scale parameters to the activity center relocation spatial scale parameters, the magnitude of the pooled‐sex, between‐year activity center relocation was similar to within‐year space use (Fig. ), suggesting that between‐year activity center shifts were not considerably larger than typical space use within years. In absolute terms, the estimated median yearly relocation distance at Gallon Jug, representative of the broadleaf sites, was 1.74 km, while the estimated median at Mountain Pine Ridge was 5.08 km. For comparison, the estimated 95% home‐range radii, computed from the detection function spatial scale parameters were 4.09 km for Gallon Jug (mean σ across sexes) and 7.76 km for Mountain Pine Ridge. Assuming study areas were perfectly circular, their mean radius was 20.75 km, and the median yearly relocation distance was 8% of the mean study area radius at the broadleaf sites and 24% of the mean study area radius at Mountain Pine Ridge. Thus, there was considerably more movement, immigration, and emigration at Mountain Pine Ridge than in the broadleaf sites. This result is confounded by the male bias at the Mountain Pine Ridge site and our inability to estimate sex‐specific activity center relocation parameters at all sites, so we cannot conclusively determine whether sex, habitat quality, or a combination of the two led to differing between‐year movements.

The broadleaf forest sites are largely homogenous in habitat, making population dynamics relatively easy to interpret; however, the population dynamics at Mountain Pine Ridge must be interpreted with respect to the surrounding landscape. The high elevation pine forest at Mountain Pine Ridge is surrounded by contiguous broadleaf forest to the south, west, and east, and more fragmented broadleaf forest to the north and northwest. The spatially explicit posterior density plots across time (Fig. ) suggested that the majority of the few individuals living on this landscape each year were located in the south portion of the pine forest, and closer to streams, or in the broadleaf forests itself. Further, given the large detection function spatial scale parameter estimate and male sex bias, this population likely consists of males that have established home ranges that overlap both pine and broadleaf forest, perhaps due to subordinate status. Then, given the large yearly relocation estimates, they may move into the pine forest in years when they cannot obtain a territory in the preferred broadleaf forest habitat and then move back into the broadleaf forest when they are more competitive. Given these dynamics, the estimated per capita recruitment rate likely includes little to no in situ recruitment within the pine forest, itself. While we estimated a declining population size through time at Mountain Pine Ridge, we do not see this as a cause for concern. The spatial dynamics illustrated in Fig.  suggests the pine forest is a marginal habitat where use varies through time, perhaps providing a refugia to subordinate males until they can obtain a territory more exclusively in broadleaf forest.

In order to correctly interpret our parameter estimates, it is important to discuss the limitations imposed by the camera‐trap, capture–recapture data. First, because juveniles with their mothers are not reliably photographed and/or identifiable, our density and survival estimates do not include this component of the population. Therefore, the density and survival estimates do not reflect a single age class, but a mixture of (primarily) adults and independent sub‐adults. Second, because individuals are included in the analysis only after they have left their mother at which point males begin to disperse, we cannot strictly estimate in situ recruitment (as in Nichols and Pollock )—juvenile males may leave their natal site before they are detected. Therefore, in order to interpret our recruitment estimates as in situ, we must assume that juvenile dispersal into and out of a site is symmetric, which may be plausible in the relatively stable broadleaf forest sites, unless they act as population sources. Then, in situ recruitment in the population source will be confounded with asymmetric dispersal out of the site and thus underestimated. Per capita recruitment should then be interpreted as the number of individuals that are born in and remain at a focal site each year plus the number that are born at another site and disperse into the focal site before they are detected via a camera trap divided by the total abundance at the focal site. Both of these events include a component of juvenile survival, although the survival period is variable due to the lack of seasonality in breeding among ocelots in Belize (Sunquist and Sunquist ). Despite these caveats, we can reliably estimate the number of new additions at each site, regardless of whether they were born there, which is the relevant quantity to estimate population growth rates.

In this study, we have provided survival, per capita recruitment, and population growth rate estimates for five sites in Belize across 4–12 yr. Our findings suggest survival is high and recruitment is low, consistent with a K‐selected species. We found evidence of a small population decline at two of four broadleaf sites; however, given the level of uncertainty, we suggest continued monitoring of these sites (Cockscomb and La Milpa) to more precisely estimate the population growth rates. Research and conservation efforts are critical to long‐term population sustainability of elusive carnivores like ocelots and other sympatric cat species, especially in areas of the world such as Belize, where deforestation, growing human populations, and infrastructure development are increasing. Estimating parameters such as population density and survival, and tracking population trends through time, are crucial to effective management of wildlife populations, particularly because these species directly influence the food web through their role as top‐down predators on smaller prey. Furthermore, we highlight the importance of estimating demographic parameters such as survival, per capita recruitment, and population growth rate as they are important determinants of population dynamics. In addition, we encourage other ecologists to collect data at larger spatial and temporal scales, to produce estimates with increased precision and decreased bias, and to better understand population dynamics of long‐lived, threatened species.

Acknowledgments

We thank the Belize Forest Department, Program for Belize, Gallon Jug Estate, Yalbac Ranch and Cattle Company, Belize Audubon Society, and Bull Run Farm (G. W. and M. Y. Headley) for permission to conduct this research. We thank the Philadelphia Zoo, Virginia Tech, Wildlife Conservation Society, and Panthera for funding over the years. We thank Richard Chandler for providing an MCMC script, which served as the foundation for the script used in our analyses. We thank Nicodemus Bol, Mauricio Aguilar, Thomas Mc Namara, former graduate students Adam Dillon, Claudia Wultsch, and Miranda Davis, and numerous volunteers and assistants who conducted field and laboratory work, from multiple universities, especially the University of Belize, The Environmental Research Institute, and Virginia Tech.