Integrating counts, telemetry, and non‐invasive DNA data to improve demographic monitoring of an endangered species

Abstract

Population monitoring can take many different forms, and monitoring elusive and endangered species frequently involves a variety of sparse data from different sources. Small populations are often hard to sample precisely and without bias, so when estimates of vital rates like survival or recruitment point to conflicting population trends, it can be hard to determine which is more correct. Furthermore, data can be extremely hard to collect on small populations and it can be helpful to find a way to use all available hard‐won data. To address these issues, we developed an integrated population model (IPM) using all available data to estimate vital rates and abundance for a case study of an endangered woodland caribou (Rangifer tarandus caribou) population. This IPM allowed us to incorporate data from juvenile recruitment surveys, telemetry‐based survival, aerial population counts and mark–resight data, and non‐invasive capture–recapture DNA data to better understand population status and trend. We estimated survival, abundance, and recruitment of four age classes of male and female caribou: young, juveniles, subadults, and adults. The four‐age class structure of the IPM allowed us to estimate recruitment from reproductive‐aged female caribou alone, even though it can be difficult to distinguish age classes—and even sexes—in the field. As part of our IPM, we developed a novel mixture model to break apart data from different age classes when age is unobservable, as it typically is from non‐invasive DNA samples. This helped us decrease bias in juvenile and adult survival estimates from scat data, which was important to our understanding of the population dynamics. Overall, our integrated model provided more precise estimates of population trends than any one method (e.g., telemetry or non‐invasive DNA) alone. This IPM provides a useful, flexible tool for biologists to monitor populations and provides a valuable example of the benefits of integrated population modeling approaches for endangered species management and recovery.

Introduction

Endangered species conservation is a growing challenge for applied ecologists and agencies as humans continue to negatively impact biodiversity (Clark et al. 2002, Schwartz 2008). Monitoring population trend and population response to management actions is key to conservation efforts (Schwartz 2008). Unfortunately, it is often difficult to make robust inferences about rare species, precisely because they are rare and sample sizes are therefore small (Schaub et al. 2007). Due to limited resources for monitoring, empirical data on many endangered species are often fragmented across time and space.

One high-profile example of an endangered species where inference is limited is the woodland caribou (Rangifer tarandus caribou). Woodland caribou are listed as endangered or threatened across much of their range in the USA and Canada (Ray et al. 2015) and are the most recent mammal to become extinct in the lower 48 states (Robbins 2018). Causes of the declines vary across space, but generally they are linked to apparent competition (Holt 1977, Hebblewhite et al. 2010, Bergerud et al. 2015, Serrouya et al. 2019). In Jasper National Park, past management practices induced apparent competition by directly altering abundance of ungulates and predators (Bradley and Neufeld 2012). As a result, caribou numbers have declined dramatically and most subpopulations of caribou in the park are extirpated or considered functionally extinct (DeCesare et al. 2011, Neufeld 2017, McFarlane et al. 2018).

Monitoring caribou population trends is a major challenge rangewide. Most woodland caribou populations occur at low density and in dense coniferous forests that make aerial surveys impractical (Courtois et al. 2003). Mountain ecotypes of woodland caribou, such as those in Jasper National Park, are somewhat of an exception because they can be surveyed above treeline during the rut (Serrouya et al. 2017), but their low densities make surveys expensive. Inclement weather and prohibitive cost can make population surveys inconsistent across time, and the amount of data collected each year tends to be small.

For species where monitoring is challenging, biologists often attempt different methods over the years to estimate population size and trend (e.g., counts, radio-collar or tag data, or citizen science data). Many other caribou populations have been monitored in the past through a simple population model, the recruitment–mortality (R-M) equation (Hatter and Bergerud 1991), which combines a marked sample of adult female caribou with aerial surveys of juvenile recruitment (DeCesare et al. 2012). More recently, researchers have begun estimating caribou population trends using non-invasive DNA from feces (Hettinga et al. 2012, McFarlane et al. 2018, 2020). Similarly, as costs of obtaining genotypes from non-invasive DNA samples have declined, DNA has become a popular monitoring tool for bears (Ursus spp., Woods et al. 1999, Kendall et al. 2009), European wildcats (Felis silvestris, Kéry et al. 2010), and dozens of other species (Mills 2012). Many of the species that are monitored through non-invasive DNA sampling are also monitored in other ways, such as through harvest records, aerial counts, or occupancy surveys.

When multiple methods are used to monitor a single population and are analyzed independently, the results from each method can be impossible to reconcile. For example, Lorenz and Raphael (2018) found decreasing densities of adult Marbled Murrelets (Brachyramphus marmoratus), stable juvenile densities, and stable juvenile:adult ratios over an 18-yr study. These results suggested high productivity, which seemed to conflict with a concurrent telemetry study showing low breeding propensity in adults (Lorenz et al. 2017). The authors proposed four different hypotheses to explain these results, demonstrating that at times, it can seem like more data actually lead to more uncertainty about the system. When different data streams conflict, the temptation is often to refer to one established method (e.g., aerial surveys or telemetry) or the more precise estimate as the gold standard (Serrouya et al. 2017). However, truth is always unknown and therefore accuracy is difficult to compare across multiple methods. This problem is widespread across species and settings, but it is especially evident with endangered or rare species. Small sample sizes result in high sampling variance and greater uncertainty, and assumptions such as independence are frequently violated (Schaub et al. 2007, Zipkin and Saunders 2018). Thus, it can be impossible to know what management decisions are appropriate and warranted when different datasets point to contrasting population trends.

Integrated population models (IPMs) can help disentangle contrasting population trends suggested by independent data sources (Saunders et al. 2019). They also can use all the available, hard-to-collect data on rare species. Integrated population models use information from multiple data sources to estimate demographic parameters simultaneously, thus creating a single comprehensible estimate of population dynamics (King et al. 2008, Schaub and Abadi 2011, Zipkin and Saunders 2018). Furthermore, IPMs use information from all the available data on a population, which usually produces more precise parameter estimates than those from independent analyses (Schaub and Abadi 2011). For example, an IPM for elk (Cervus elaphus) showed that combining telemetry with aerial survey data reduced the uncertainty in demographic parameter estimates by up to 30% (Eacker et al. 2017). Integrated population models can then be used to assess tradeoffs between the cost and effort associated with each dataset against the precision increase it contributes, which is particularly important in monitoring programs for endangered species (Ahrestani et al. 2017, Sanderlin et al. 2019).

An IPM is a hierarchical model comprised of a biological process model and a set of observation models (Schaub and Abadi 2011). First, the biological process model is a demographic model that defines how the biological parameters (e.g., survival, recruitment, and abundance) connect to each other. For example, $N_t=N_{t-1}(S_t+R_t)$ is a simple biological process model, in which N_t, abundance at time t is a function of abundance at the previous time step N_t−1, survival S_t, and recruitment rate R_t. The process model can incorporate additional complexity by adding age classes, sexes, age of first reproduction, or time-varying parameters. The process model should capture the relevant complexity of the system without adding too many parameters for the available data. The second part of an IPM is the set of observation models. Observation models link an incoming data source to a biological parameter (Schaub and Abadi 2011). For example, a known-fate survival observation model takes the observed states (live or dead) for a collection of individuals tracked over time to estimate a period-specific survival for the population on average (S). The survival parameter S was also defined in the biological process model, thus establishing the link between the process model and observation model. An observation model is built for each available data source (e.g., in addition to the known-fate observation model, a second observation model could map observed young:adult ratios to another biological parameter of interest, recruitment, Harris et al. 2008). The two-tiered structure of the IPM (biological process model and set of observation models) allows all biological parameters to be estimated simultaneously and within the bounds of the biological process model. The structure integrates data streams to estimate one single population trend and results in parameter estimates that are more precise than those from a single data set alone. Because the biological process model can be very complex, IPMs can also estimate sex- and age-specific demographic rates and immigration/emigration rates, even when such data were not collected (Schaub et al. 2010), although it is important to consider the assumptions of IPMs when interpreting latent parameters (Riecke et al. 2019).

Although managers often monitor populations with DNA and other methods concurrently, approaches for integrating information between these methods are largely unexplored. In particular, for species with complicated life histories (e.g., late age of first reproduction) that require complex stage-structured models, it is critical to separate data for each life stage. Although genotypes provide a unique ID, the age class of the individual is usually unknown (Furnas et al. 2018). While promising techniques are emerging to estimate age from fecal samples, their use is not yet widespread and further statistical development is needed to account for potential error (Haussmann and Vleck 2002, Flasko et al. 2017, McFarlane et al. 2018). We developed a novel way to incorporate non-invasive DNA data from individuals of unknown age into a complex stage-specific IPM. This advance allowed us to improve estimates of juvenile recruitment and adult survival, and these key vital rates are often especially important in declining or endangered populations.

Managers of caribou populations are facing a crisis known to many working with endangered species; decisions need to be made quickly, yet they need to be informed by accurate and precise information. When debating drastic actions like implementing conservation breeding or translocations, precise and accurate results are critical but extremely difficult to obtain. It is important to have a way to integrate each hard-earned data point into a flexible, comprehensive population monitoring tool that produces fast results and allows users to project recovery scenarios based on real-world data. When assessing potential management actions and modeling future scenarios, wildlife managers can make better decisions when debate over conflicting datasets is minimized and tools to integrate the data are readily available. To be useful, this model must be fast and user-friendly; managers cannot wait months for an MCMC chain to converge. Here, our goal was to develop a fast, flexible IPM to estimate caribou population caribou trends, using available data from aerial population estimates, juvenile recruitment surveys, radio-collared adult females, and non-invasive DNA mark–recapture surveys (McFarlane et al. 2018). Our approach produced a picture of caribou population dynamics and also offers general insight for integrating non-invasive DNA data with other, more conventional data such as aerial surveys or telemetry.

Methods Study area

The study area is located in the southern half of Jasper National Park (JNP), Alberta, Canada (Fig. 1; 52°8′–53°21′ N, 116°48′–119°27′ W). The mountainous topography ranges from 1000 to 3500 m in elevation. The climate is characterized by long, cold winters and short summers, with most precipitation occurring in spring. The landscape is ecologically classified into the montane, subalpine, and alpine ecoregions (Holland and Coen 1983). In addition to the main predator of caribou, gray wolves (Canis lupus, Hebblewhite et al. 2007), the predator community also includes mountain lions (Puma concolor), grizzly bears (Ursus arctos), black bears (Ursus americanus), and wolverine (Gulo gulo).

Fig. 1. Caribou subpopulations in Banff and Jasper National Parks. The Banff and Maligne subpopulations have been extirpated, and the Brazeau subpopulation is considered functionally extinct. Source DeCesare et al. (2011).

Caribou within the Rocky Mountain national parks are federally classified as the threatened Southern Mountain population (Species at Risk Act 2002, Ray et al. 2015), although they are recommended by COSEWIC to be re-designated as Central Mountain caribou, and up-listed to endangered status (Committee on the Status of Endangered Wildlife in Canada [COSEWIC] 2014). The present Southern Mountain caribou recovery strategy (Environment Canada 2014) applies to the caribou in Jasper National Park and identifies the Banff, Brazeau, Maligne, and Tonquin subpopulations as the Jasper/Banff local population unit (LPU). The Banff and Maligne subpopulations are considered extirpated, and the Brazeau subpopulation is considered functionally extirpated (counts of fewer than 10 reproductive females), so here we focus on the Tonquin subpopulation in southern Jasper (Fig. 1).

Study species: caribou

Caribou are unique in the cervid family because females have antlers, making it challenging to discriminate sexes during aerial surveys. In comparison to other cervids, caribou also have relatively slow reproduction. Caribou typically do not reproduce until their third year of life (Fancy et al. 1994, Adams and Dale 1998) and they have only one young per reproductive event. Twins occur at a rate of less than one in 1000 reproductive events (Hummel and Ray 2008). As with other ungulates, adult female survival is high and relatively constant, whereas juvenile survival is variable and often lower than adult female survival (Gaillard et al. 2000). While there has been only one empirical study of radio-collared woodland caribou juvenile survival to 1 yr (Gustine et al. 2006), dozens of studies using age recruitment data (the number of young per 100 females in late winter) have confirmed that juveniles have lower survival than adults (DeCesare et al. 2012, Hervieux et al. 2013). Due to these life-history characteristics, it is critical for population models to capture the differences in fecundity and survival rates between juveniles and adults.

Biological process model

The biological process model in an IPM defines how the demographic parameters of survival, abundance, and reproduction relate to each other between years. Due to limited movement between caribou subpopulations (van Oort et al. 2011, McFarlane et al. 2018; L. Neufeld, unpublished data), we did not model immigration or emigration in this population. To account for differing survival between ages and sexes and the delayed age of first reproduction, we structured the caribou biological process model into four age classes for each sex. As is typical with many ungulates, our data were all collected during the autumn, so we defined the model year to occur between October 1 and September 30. On October 1 of year t, individuals fell into the following age classes: young (0.5 yr), juvenile (1.5 yr), subadult (2.5 yr), and adult (3.5+ yr). We defined survival for each age class as the period from t − 1 to t. For example, juvenile survival referred to individuals during the year they were from 0.5 to 1.5 yr of age.

We modeled the abundance (N) of juveniles, subadults, and adults in year t as a function of abundance the previous year and survival (S). As an example, the number of juveniles (J) of sex s in year t was modeled as[Image Omitted. See PDF]

In our implementation, we used a Normal approximation of the binomial distribution, which allowed for faster run times and better MCMC chain mixing (Brooks et al. 2004). This approximation took the form:[Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF]

Subadult and adult abundance were modeled similarly to juvenile abundance, using abundance of the appropriate age class in time t − 1. Adult abundance depended on both subadult and adult abundance in year t − 1.

To appropriately fit our data and capture biological reality, we defined both an average recruitment rate for all age classes (RA_t: calves per female 1.5 yr and older) and a true recruitment rate for reproductive-aged caribou (R_t: calves per female 2.5 yr and older). It is nearly impossible to differentiate female juveniles, subadults, and adults reliably by sight in the field. Therefore, in the field we collected ratios of young to female caribou aged 1.5 yr and older. These observations were used to estimate RA_t, the average contribution of juvenile (J), subadult (S), and adult (A) female caribou to recruitment. We matched the biological process model to fit this interpretation. A major advantage of the IPM is that it estimates abundance of each age class separately. Therefore, we also were able to derive true recruitment, R_t, as the ratio of young to reproductive-aged females (2.5+ yr).

We estimated abundance of young (Y) caribou from the abundance of female (F) caribou and the average recruitment rate RA_t. Our observations occurred in the autumn, so RA_t encompassed pregnancy rate, mother survival, and neonatal survival to six months of age. To model abundance of young caribou, we used a Normal approximation of the Poisson distribution for animals of sex s in year t:[Image Omitted. See PDF]

The Poisson distribution allows for the possibility (however rare) of twinning, which the Binomial distribution would not allow. In the absence of evidence suggesting that neonatal sex ratios are skewed in caribou, we divided the recruitment term in half to calculate the number of male and female young.

Observation models

Observation models in an IPM connect each available data source with the biological parameters defined in the process model. In our model, these parameters were survival, recruitment, and abundance for each of the four age classes and each sex. We collected up to five data types for the Tonquin subpopulation annually between 2003 and 2019: (1) Lincoln-Petersen estimates from aerial surveys, (2) age ratios from aerial surveys, (3) telemetry-based survival, (4) non-invasive scat DNA-based capture–recapture data, and (5) minimum known alive counts (Fig. 2). Not all data types were available for every year, which we note in the description of each data type. We provide a brief summary here of data collection methods and refer readers to the details provided by Hebblewhite et al. (2007), DeCesare et al. (2011), Bradley and Neufeld (2012), McFarlane et al. (2018), and Neufeld and Bisaillon (2017).

Fig. 2. Visual representation of available data (green diamonds) for the caribou integrated population model in Jasper National Park. The demographic parameters (red rectangles) and uncertainty (blue rectangles) are estimated together in the integrated population model from the available data.

We implemented a two-stage approach to fitting the IPM to our data (Besbeas et al. 2002). We analyzed each of the five datasets independently, which produced annual estimates and their standard errors for each dataset. Put together, the estimates and standard errors represented the distributions of possible values we could have collected for each dataset in each year. To incorporate these into the IPM, we created an observation model for each dataset that drew values from the dataset's distribution. This technique allowed very fast convergence and chain mixing, which can be very difficult to obtain with IPMs. This technique propagates error from each dataset and provides similar inference to incorporating the raw data into the IPM likelihood. Two of our datasets were straight-forward analyses: Lincoln-Petersen and age ratios. We used a Normal distribution in the IPM to draw values from the estimated data distribution. We used a similar approach on two of our other data types—telemetry and non-invasive scat DNA—but we go into additional detail on their analysis due to additional complexity. Our final data type, minimum counts, did not have standard errors, so we built an observation model in the IPM that incorporated the raw data directly into the likelihood.

Lincoln-Petersen

We performed aerial surveys every fall (late September, early October) from 2003 to 2019. From 2003 to 2014 the number of marked (collared) and unmarked individuals seen were used for mark–resight Lincoln-Petersen estimates of abundance (Chapman 1951). We were not able to estimate abundance using this method in 2005 and 2012 due to a lack of resighted collars and inclement weather, respectively. Misidentification can be common in caribou ratio data because both males and females can have antlers at the time of the survey. To address this, we classified animals in the field by using multiple observers and image-stabilizing binoculars, as well as by landing and using a spotting scope on larger groups. We also had experienced photographers take high-quality digital photos of each individual for later confirmation of sex, which greatly increased our confidence that the adult female category did not mistakenly include male caribou.

Our IPM observation model for the Lincoln-Petersen abundance estimates was a Normal distribution. This allowed for process variation around the point estimate. We drew each Lincoln-Petersen estimate m_t from the total abundance N_t and the Lincoln-Petersen standard error SE(m_t), in the form[Image Omitted. See PDF]

Age ratios

We observed ratios of young caribou to females aged 1.5 yr and older through our aerial counts. Similar to the Lincoln-Petersen observation model, we used a Normal distribution to estimate the contribution of juvenile, subadult, and adult female caribou to average recruitment (RA) using the observed age ratios r and their standard error:[Image Omitted. See PDF]

Known fate survival model

Thirty-four adult female caribou were captured and equipped with radio-collars from 2003 to 2011. These collars were monitored for alive/dead status from the air every 1–3 months until 2017, a monitoring design that provides an unbiased estimate of survival for adult female caribou under these conditions (DeCesare et al. 2016).

Before incorporating these data into the IPM, we estimated adult female survival with a known fate model. We estimated survival φ in month m from the monthly fates for each animal (y_m):[Image Omitted. See PDF]

We included a month-of-year random effect on survival to account for cyclical differences in survival throughout the year. We derived annual survival φ_m by multiplying the 12 months' survival together.

To bring this data source into the IPM, we fit a Normal distribution from the mean estimated survival and its standard error to inform the survival parameter S for adult females:[Image Omitted. See PDF]

DNA sampling and analysis

Fecal samples were collected from the Tonquin in JNP over eleven consecutive winters from 2007 to 2018. We collected pellets in a robust design, with two pellet surveys each year, spaced at least 10 d apart (Pollock 1982, Kendall et al. 1997). This resulted in 22 total sampling occasions. Individuals in each occasion were identified using 10 polymorphic microsatellite loci and a primer for sex identification. Previous studies in JNP have used non-invasively collected DNA from aerial surveys to estimate population size, trend, recruitment, survival, and individual fitness (McFarlane et al. 2018). We refer readers to Hettinga et al. (2012) for details of DNA sampling design, DNA extraction, genotyping, and construction of individual capture–recapture histories from non-invasive samples.

Outside the framework of the IPM, we developed a modified robust design model (Pollock 1982, Kendall et al. 1997) to create estimates and standard errors of survival and abundance from the non-invasive DNA samples. Although the scat samples were identified to sex, they came from animals of unknown age. While new methods offer promise to associate age with pellet samples (Flasko et al. 2017), results are preliminary and were not available for all the genotypes or years in our capture histories. This problem is common in non-invasive sampling (Flasko et al. 2017). Although subadult and adult survival is often assumed to be equal (DeCesare et al. 2012), juvenile survival is typically lower in most ungulates (Gaillard et al. 2000). Therefore, ignoring age in a robust design model may produce biased estimates of both juvenile and subadult/adult survival.

To separate juvenile survival from subadult/adult survival, we developed a model for the true, unobservable age of an individual. This latent parameter was informed by subsequent recaptures of an individual. We modeled the true age (g) of individual i in the year of its first capture by[Image Omitted. See PDF]where π_t was the probability that individual i was an adult on its first capture occasion. Therefore, we let g = 1 if the animal was an adult and g = 0 if the animal was a young. For the first encounter of an individual, we defined the survival probability in the robust design to be[Image Omitted. See PDF]

Whenever that individual's genotype was detected in subsequent years, we knew the animal must have been at least 1.5 yr old. Assuming equal subadult and adult survival, we assigned the individual the adult survival S_A in every year following its first capture. We allowed π_t to vary over time to avoid the effects of mark saturation (i.e., probability that all adults have been previously captured at least once).

We ran two chains for 50,000 iterations after 10,000 iterations of burn-in using JAGS version 4.2.0 (Plummer 2003) and R version 3.5.3 (R Core Team 2019). Convergence was assessed through visual inspection and R-hat values <1.1. We incorporated the resulting estimates of survival (f_t,s,a) and abundance (d_t,s,a) and their standard errors into the IPM using a Normal distributions, similar to the other observation models:[Image Omitted. See PDF][Image Omitted. See PDF]

Minimum counts

Minimum counts were made using all available data on uniquely identifiable individuals. These typically came directly from aerial surveys, but in some cases we collected scat DNA from additional collared individuals who were not seen in the aerial survey but whose genotypes were known. Minimum counts (c) served as lower bounds for our IPM estimates of abundance. The observation model, therefore, was a negative binomial distribution:[Image Omitted. See PDF][Image Omitted. See PDF]

The negative binomial distribution with abundance as a parameter specifies a lower bound below which abundance cannot go.

Model implementation

The joint likelihood for the observation models was the product of the likelihoods of each component model:[Image Omitted. See PDF]

We used informative priors (standard deviation 0.1) to capitalize on our knowledge of caribou population dynamics in this system based on previous studies (Hebblewhite et al. 2007, DeCesare et al. 2011, Bradley and Neufeld 2012). For recruitment, the informative prior centered around the mean observed age ratio (0.41). For juvenile survival, our informative prior centered on 0.65, and for subadult and adult survival it centered on 0.85, based on observed survival for other caribou populations (Gustine et al. 2006, Hebblewhite et al. 2007, DeCesare et al. 2011, Hervieux et al. 2013). The prior on the first-year population size for each sex and age was derived from the first estimated population size from the Lincoln-Petersen estimates. We ran three chains for 100,000 iterations after 50,000 iterations of burn-in using JAGS version 4.2.0 (Plummer 2003) and R version 3.5.3 (R Core Team 2019). We assessed convergence by visual inspection and R-hat values <1.1. We allowed recruitment, juvenile survival, and adult survival to vary over time by adding a random effect for year. We chose to let these parameters vary based on our understanding of the underlying biology of caribou populations at small numbers. Variability in adult female survival is rarely observed in long-lived ungulates, which is consistent with the fact that these populations are highly sensitive to changes in adult female survival (Gaillard et al. 2000). Thus, as abundance approaches zero, the death of any one female has the potential to greatly influence population trajectory. We provide annotated JAGS code for our integrated population model for caribou, including the novel integration of non-invasive CMR data (Data S1).

Finally, we tested goodness of fit in a two-staged approach. First, we tested goodness of fit of each dataset to its observation model independently, if there was an available test. Only our robust design scat data had an available test. We calculated its chi-square goodness of fit using Test 2 and Test 3 (Burnham et al. 1987) in program RELEASE (RMark 2.2.7, Laake 2013). Then, we calculated an overall goodness of fit of the IPM abundance estimate to the observed scat abundance dataset by calculating a Bayesian P-value (Gelman et al. 1996, Besbeas and Morgan 2014).

Results

The IPM structure allowed us to estimate survival and abundance of each sex and age class, as well as recruitment and the derived population growth rate, λ (Fig. 3). We derived λ from the estimates of abundance in years t and t − 1. Sample sizes for all data types are reported in Appendix S1.

Fig. 3. Integrated population model demographic estimates for Tonquin caribou, including (a) female abundance (aged 2.5+ yr), (b) adult female survival (aged 3.5+ yr), (c) recruitment, and (d) population growth rate (λ). The dashed line at λ = 1 represents stable population (no growth, no decrease). In each panel, the heavy line is the median estimate, and the shaded area represents the 95% Bayesian credible interval.

The median estimated adult female survival was 0.703, with a 95% Bayesian credible interval (CRI) of 0.639–0.77 (Fig. 3). Median subadult female survival was 0.75 (95% CRI 0.73–0.77; Appendix S2). Juvenile female survival was highly variable between years, but the median was 0.699 (95% CRI 0.21–0.98, Appendix S2). Over the first few years of monitoring, sample sizes were small on many of our data sources, which is reflected in the wide 95% Bayesian credible intervals on the IPM estimate. When we began monitoring in 2003, the population was already depressed in comparison with historical lows, and it showed a continual decline in abundance from 2008 forward. The population was estimated to be about at the functional extinction threshold of 10 females aged 2.5 yr and older in 2017 onward.

The four age classes of the IPM biological process model allowed us to estimate abundance of each age class, regardless of whether we could reliably distinguish them in the field. Even though recruitment data were collected as observed ratios of young to female caribou aged 1.5 yr and older, we were able to estimate true recruitment, R_t, as the ratio of young to breeding-aged females (2.5+ yr) in a given year (Figs. 3, 4). This is a marked difference from the raw field ratios where juvenile, subadult, and adult female caribou cannot be distinguished, which only allows estimation of the average contribution of females age 1.5 yr and older to recruitment (Fig. 4).

Fig. 4. Comparison of integrated population model (IPM) estimates of recruitment (gray line) with 95% Bayesian credible intervals (shaded region) against observed age ratios (red points, with 95% confidence intervals). Because caribou age classes are difficult to tell apart in the field, the red points represent ratios of young to female caribou aged 1.5 yr and older. In contrast, the IPM estimate was derived from the ratio of young caribou to females of reproductive age only. Due to inclement weather, we were not able to obtain an age ratio in 2012.

The integrated nature of the IPM allowed us to investigate the relative importance of each dataset on the final estimate. The IPM takes information from all available data and weights them by their standard error, which results in the more precise data more heavily influencing the final estimates. For instance, the non-invasive DNA abundance estimates had much smaller standard errors than the Lincoln-Petersen abundance estimates and therefore appeared to drive the IPM abundance estimates (Fig. 5).

Fig. 5. Total abundance estimate from the integrated population model (IPM; gray line) with 95% Bayesian credible intervals (shaded region), plotted with the Lincoln-Petersen abundance estimates and their 95% confidence intervals (red) and the non-invasive DNA abundance estimates with their 95% credible intervals (blue). The IPM closely tracks the estimates from the non-invasive DNA model because the standard errors are much smaller than those of the Lincoln-Petersen estimates.

To test goodness of fit, we first conducted a chi-square goodness of fit test for our robust design scat data. We found no evidence of lack of fit for this observation model (P = 0.89). Then, we tested the IPM abundance estimate's goodness of fit to our observed scat abundance, and we found poor overall fit (P = 0.02). We investigated yearly goodness of fit and found lack of fit in 2009, when scat abundance confidence interval fell outside the IPM estimated abundance (Table 1). We found improved model fit when we removed 2009 from the P-value calculation (P = 0.07), and even better fit after removing both 2008 and 2009 (P = 0.19). Although the cutoff of P < 0.5 is arbitrary, this generally demonstrated that 2009, and to a lesser extent, 2008, drove the overall lack of fit.

Table 1 Bayesian P values showing goodness of fit for the integrated population model by year.

Year	P
2007	0.424567
2008	0.11069
2009	0.02098
2010	0.240443
2011	0.420463
2012	0.48937
2013	0.496833
2014	0.1143
2015	0.50006
2016	0.49174
2017	0.49231
2018	0.430553

Note:

Results demonstrate that the overall fit of the model is strongly influenced by a lack of fit in 2009 (P < 0.05).

Discussion

Our IPM showed a steady decline in caribou abundance in the Tonquin from approximately 101 individuals in 2008 to 38 in 2019. These results are consistent with the high number of noted extirpations in the last decade at the southern edge of caribou range in western North America (Hebblewhite et al. 2009, Robbins 2018). Combining all available data on the Tonquin caribou subpopulation into an integrated population model allowed us to estimate sex- and age-specific survival, abundance, and recruitment concurrently for this imperiled population. The integrated nature of the model painted a single coherent picture of the status and trend of the population. We found overall low adult female survival along with highly variable juvenile survival and recruitment. Adult female survival is a highly elastic demographic rate in long-lived ungulate species, and low values can lead to dramatic population declines (Gaillard et al. 2000).

A previous study on this population used Lefkovich matrix modeling to project consequences of different translocation scenarios, and the authors concluded that the future viability of the Tonquin population was high without additional recovery actions (DeCesare et al. 2011). Our results, in a different analysis framework and with nearly a decade more data, instead suggest declining viability of the Tonquin herd, commensurate with especially low adult female survival and recruitment rates. Our work and previous studies indicate that urgent recovery action is needed to reverse the declines the Tonquin population has experienced (Hebblewhite et al. 2007, DeCesare et al. 2011, McFarlane et al. 2018). The number of remaining breeding females is likely too low to support recovery without augmentation of their numbers (Gilpin and Soule 1986, Hervieux et al. 2013). Because the population decline we recorded was likely driven by very low adult female survival, population augmentation will have to be accompanied by improvement of the conditions that led to poor survival during our study (Hebblewhite et al. 2007). We were limited in our ability to fully explore the drivers of and variation in adult female survival, and we did not have the resolution to detect whether female survival may already be improving due to management interventions. However, the flattening rate of decline in recent years may indicate that conditions have improved, which we will further investigate in future work.

Our model builds on previous demographic models of woodland caribou that have used the simple R-M equation with only information about adult female survival and juvenile recruitment (DeCesare et al. 2011). While many populations of woodland caribou and other ungulates are effectively monitored via this simple approach, there is an increasing effort to collect concurrent data using alternative methods. Serrouya et al. (2017) compared population trend estimates from aerial counts to demographic modeling via the R-M equation yet no studies have formally combined different methodological approaches into an IPM to coherently estimate trend in caribou. The IPM approach allowed us to use multiple data sources for a small population of endangered animals, and it helped us navigate potentially conflicting datasets. A significant benefit of the IPM over other age-structured population models that have been used to assess viability of caribou (DeCesare et al. 2011) is that IPMs do not assume a stable age distribution. The assumption of stable age distribution is usually required for most applications of Leslie matrices, or, stage-age-based Lefkovich models, for example, those used by DeCesare et al. (2011) to previously predict caribou viability in Jasper. Integrated population models, however, make no such restrictive assumptions of age structure, and instead allow the data to inform time-varying age structure. This was particularly helpful in our analysis for the Tonquin population because it allowed us to assess recruitment and survival as time-varying parameters.

Our approach also allowed us to infer age information from non-invasive DNA when combined with age ratio and other data sources. A common problem of non-invasive DNA surveys is that the age of the individual animal that deposited a fecal sample is typically unknown. For example, non-invasive monitoring of gray wolves (C. lupus) in France and Spain (Cubaynes et al. 2010, López-Bao et al. 2018) and grizzly bears (U. arctos) in Alberta (Boulanger et al. 2004) noted substantial individual capture heterogeneity that may have been due to unknown age structure of non-invasively collected DNA samples. In contrast, our IPM used information collected via conventional surveys to help estimate juvenile survival rates in an integrated fashion. Without this approach, our estimates of adult survival would have been biased because of the unknown fraction of adults, subadults, and juveniles in the sample. Because of the strong differences in adult and juvenile survival for most species' life histories, failing to account for such heterogeneity would had led to biased estimates of demography (Boulanger et al. 2004). The general structure of our IPM also provides a strong framework for building in age information that can be determined from DNA samples, using techniques like measuring telomere length or hormones (Haussmann and Vleck 2002, Flasko et al. 2017, McFarlane et al. 2018). The age information gained from these techniques may further improve the precision of the juvenile and adult survival estimates in our IPM. Because these techniques are so new, the information is often unavailable for all years or older genotype samples, as was the case in our study. In future analyses, we could develop an additional observation model to incorporate this information and account for the sources of uncertainty.

An important part of our approach to understand the population trajectory of the Tonquin caribou was the use of informative priors. Informative priors in the IPM formalize a way to build on established knowledge (Nowak 2015). We were able to build a foundation based on the ecological principles common to long-lived, slowly reproducing species like caribou and other large ungulates (Gaillard et al. 1998). Because we were working with limited data due to the small population size of our study population, we worked with caribou experts to build the most ecologically meaningful priors for this system, based on years of knowledge built on caribou life history and parameters collected for similar ungulate species (Gaillard et al. 2000, DeCesare et al. 2011, Hervieux et al. 2013). Our priors helped guide the model when there were no data, which was important for informing parameters that are extremely difficult to measure in the field, like juvenile survival. The downside to using informative priors in the absence of data is that we can never be sure whether they reflect reality. It is a truism that the more data that are available, the fewer assumptions there are necessary to make. The inverse is also true; it can be necessary to add assumptions when few data are available. Yet in our case, our data clearly supported that adult female survival over the entire study period was indeed much lower (0.70) than previously published studies that we used for our prior (0.85), and conversely, estimated recruitment was very similar to its prior. Endangered populations pose a major challenge to traditional statistics—it is often simply impossible to collect enough data to produce a precise and informative estimate. Managers in these situations still need defensible and meaningful results; by building on knowledge based on years of literature on caribou, we produced an estimator that was biologically meaningful and useful in a management context. The IPM also has the flexibility to allow the data to direct the parameter estimate away from the prior when there are sufficient data that conflict with the prior.

All IPMs make various assumptions related to both the IPM structure and the incoming data sources (Riecke et al. 2019). As with all population models, we also made certain assumptions that may not exactly reflect reality (DeCesare et al. 2016). For example, we assumed that female caribou do not reproduce until their third year of life, while in reality some small portion of younger females may reproduce, or, alternately, delay reproduction until year 4. There are always logistical and financial constraints that prevent us from collecting data on every possible parameter, which the model should reflect. When defining the process model, we kept the concept of parsimony in mind when compromising between model complexity and the amount of available data to inform extra parameters. On the data collection side, each incoming dataset also brings its own additional assumptions to an IPM. For example, assumptions in this paper include the assumption of random capture of radio-collared animals, randomized non-invasive surveys, collars not influencing the fates or survival of collared animals, and fate of radio-collared animals being known with certainty (DeCesare et al. 2016). These assumptions are based on the sampling designs of the input data, and if violated they cannot be fixed by the IPM (Riecke et al. 2019). It is important to endeavor to meet as many assumptions as possible, although this becomes increasingly difficult as population sizes get smaller. Combining multiple datasets into an IPM introduces another set of assumptions: that all incoming datasets are independent and sample the same population (Riecke et al. 2019). When working with small and endangered populations, it is especially hard to meet these assumptions. When working with data collected in the past, it becomes even harder to correct data collection issues. Regardless, the most likely effect of violating the assumption of independent datasets is overinflated precision of the resulting estimate. Thus, we caution that our estimates may be artificially precise, and that the true credible intervals may be wider than those reported here, and that this caution will apply to IPMs developed for almost all small and/or endangered populations.

All models must strike a balance between complexity and the amount of data available to inform them. This is especially true given the challenge of Bayesian model selection (Hooten and Hobbs 2015). When developing our IPM, we explicitly tested different model structures such as fixed annual survival and recruitment rates, but ultimately selected our model structure using a combination of expert knowledge of our system and applicability to management decisions (Gelman 1996). Goodness of fit for IPMs is an area of active research, and no standard method has emerged to test it (Schaub and Abadi 2011, sensu Hooten and Hobbs 2015, Conn et al. 2018). We followed a two-step approach to goodness of fit, by first testing fit for each observation model that had a test and then calculating a Bayesian P-value for overall model fit. In the first step, we only had one observation model that had a goodness of fit test, and we found no evidence of lack of fit. In the second step, we found overall support for fit of our IPM after removing two years of scat observations. The two years in question demonstrated a motivating example for creating the IPM: When faced with two datasets that indicate different trends, which one should managers believe? From 2008 to 2009, the aerial survey indicated a dramatic decrease in abundance, whereas non-invasive DNA estimates demonstrated population stability over the same period (Fig. 5). The IPM appears to have reconciled the differences between the two incompatible datasets, resulting in the IPM being drawn away from the scat dataset, regardless of its small standard error. The IPM integrated all available data—not just the abundance datasets—weighted by their uncertainty (due to multiple factors, including sample size), and the integrated estimate shows that population stability was more likely in the period 2008–2009. Integrated population models, like all shrinkage estimators, should constrain estimates and reduce noise across individual datasets, which is the behavior we saw in this instance (Burnham and White 2002). Our integrated estimate demonstrated that a manager may not always want to pick a single dataset to make inference, because the trends in the dataset with the smaller confidence intervals may not agree with other vital rates (i.e., reproduction and survival).

We formulated our IPM to be used in a decision-making framework to evaluate the assumptions and efficacy of different recovery scenarios and help managers make challenging decisions. Nowak et al. (2018) demonstrated the application of a similar IPM for harvest management of mountain lions (Felis concolor) in Montana in a decision-making framework. As in that example, we used a two-stage observation model approach to increase speed of convergence and chain mixing, thus allowing managers to obtain results in a matter of seconds, not days or months. This was important toward our goal of creating a tool that would allow managers the ability to update models quickly and compare multiple future scenarios. In our IPM, we included the ability to assess population augmentation scenarios by adding individuals over time via translocation (DeCesare et al. 2011), conservation breeding (Bisaillon and Neufeld 2017), or maternity penning (Smith and Pittaway 2011), with the ability to account for post-translocation survival depression (DeCesare et al. 2011). Although comparison of future potential scenarios was beyond the scope of this paper, it was a key driver in building our fast, flexible IPM that could be used for assessing potential future scenarios and making management decisions. Beyond evaluating future scenarios, IPMs can also be used to test hypotheses about the mechanisms of population declines through analysis of covariates' effects on demographic rates. Previous studies in our system clearly identified unsustainable rates of predation, largely from wolves (Hebblewhite et al. 2007, Bradley and Neufeld 2012), which was facilitated by roads and trails (Whittington et al. 2011), as the proximate drivers of caribou declines. Jasper National Park has implemented mitigation measures to address these factors over the course of the study period, but demographic stochasticity and Allee effects due to small population sizes (DeCesare et al. 2011) made it challenging to explicitly test these complex covariates in our analysis. However, the inherent ability to test covariates and evaluate recovery actions make this IPM potentially useful for a wide variety of other rare or endangered species.

The integrated nature of this analysis means that we created a dynamic model, which can easily be updated with additional data or covariates in future years. We believe a distinct advantage of our model formulation is that incoming data can help inform previous estimates in the time series, especially in years when data were sparse. The additional data help strengthen overall inference from the model. However, we recognize that it can present a challenge to managers communicating with the public when point estimates in the past can change slightly with the addition of new data. We suggest that a good practice when interpreting IPM estimates is to draw inference from the posterior distribution as a whole, rather than the exact point estimates. Embracing structural uncertainty is a good general practice in ecology and is essential when working with endangered populations.

The IPM developed here formalized a way to use all available data on an endangered population to estimate sex- and age-specific vital rates. Our IPM provides a flexible and customizable framework (and JAGS code) for use in other populations of caribou or other endangered species where monitoring includes multiple types of data collection. It was founded on an ecological process model specific to the species and allowed estimation of observation error inherent in data collection. When we had multiple datasets to inform a single parameter, such as the non-invasive DNA data and collar mark–recapture data to inform abundance, the IPM allowed the data with smaller standard errors to guide the estimates more strongly. Our approach will have an important impact on caribou management in Jasper National Park and can serve as the basis for a flexible tool useful in other caribou populations and the broader field of monitoring endangered species.

Acknowledgments

Funding was provided by Parks Canada, Jasper National Park, University of Montana, Alberta Environment and Parks, and NASA's Arctic Boreal Vulnerability Experiment program (ABoVE) to M.H. (NNX15AW71A). The authors have no conflicts of interest to declare.