An assessment of red wolf (Canis rufus) population persistence using classical and individual-based modeling approaches by Peter John Mahoney A thesis submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Master of Science Auburn, Alabama May 9, 2010 Keywords: Canis, wolf, coyote, demographics, hybridization, population viability analysis, PVA Copyright 2011 by Peter J. Mahoney Approved by Todd Steury, Chair, Assistant Professor of Forestry and Wildlife Sciences James Grand, Professor of Forestry and Wildlife Sciences Michael Wooten, Professor of Biological Sciences Abstract Population viability analysis (PVA) has become the quintessential analytical tool for conservation biologists. However, as interest in PVA has grown, many notable concerns have been raised about the predictive accuracy of PVA models. Model uncertainty, the analytical limitations of sparse data, and error propagation are a few of the many limitations discussed in the literature. Unfortunately, PVA models are rarely re-evaluated or tested relative to future trends within populations, limiting our understanding of how robust PVA can be to many of the uncertainties raised by the conservation community. Here, we explore how model uncertainties and model complexity influence predictions of population persistence with twenty one years of red wolf (Canis rufus) demographic data in an effort to better understand PVA, as well as guide more informed red wolf management in the future, using both classical and individual-based approaches. ii Acknowledgments I would like to thank my advisor, Dr. Todd Steury, and my committee, Dr. James Grand and Dr. Michael Wooten, for tolerating my not always so clear discussion of a complex topic. Thank you for your advice, support, and direction during my time here at Auburn. I would also like to thank the US Fish and Wildlife Service, Dr. Dennis Murray and the Red Wolf Recovery Team, and the USFWS red wolf recovery field personnel for providing the data, resources, and incentive for this population modeling project. Last but not least, to my wife Heather and daughter Gwynn, as well as to the Mahoney and Trew clans, for without them this would never have been possible. iii Table of Contents Abstract.........................................................................................................................................ii Acknowledgments........................................................................................................................iii List of Tables................................................................................................................................vi List of Figures.............................................................................................................................vii List of Abbreviations..................................................................................................................viii Chapter 1 : Effects of model structure and time series length on population viability predictions in the red wolf (Canis rufus) ................................................................................1 Introduction .....................................................................................................................1 Methods ...........................................................................................................................4 Kelly et al. model description ....................................................................................4 Current data and vital rate estimation ........................................................................5 Current model description and process .....................................................................6 Assessing model performance .................................................................................10 Results ...........................................................................................................................10 Model structure ........................................................................................................10 Length of data collection period ..............................................................................12 Perturbation analysis ...............................................................................................13 Discussion ......................................................................................................................14 iv Model structure ........................................................................................................14 Length of data collection period ..............................................................................15 Chapter 2 : The role of behavior and management in predicting red wolf (Canis rufus) persistence potential: an individual-based modeling approach ............................................26 Introduction ...................................................................................................................26 The Model .....................................................................................................................29 Purpose of model .....................................................................................................30 Entities, state variables, and spatial scales ..............................................................30 Process overview and scheduling ............................................................................32 Design concepts .......................................................................................................32 Sub-models ..............................................................................................................34 Simulations & Results ...................................................................................................38 Model robustness .....................................................................................................38 Management effort and extinction risk ....................................................................43 Sterilization v.s. euthanasia .....................................................................................45 Hybridization rates ..................................................................................................47 Elasticity analysis .....................................................................................................48 Discussion ......................................................................................................................49 References .................................................................................................................................67 v List of Tables Table 1.1 Comparison of model structure...................................................................................20 Table 1.2 Parameter values for each model.................................................................................21 Table 1.3 Performance statistics for each model.........................................................................22 Table 1.4 Perturbation analysis results........................................................................................23 Table 2.1 Model parameters and sub-models..............................................................................54 Table 2.2 Robust analysis results................................................................................................55 Table 2.3 Canid proportions inside and outside of core assessment area...................................56 vi List of Figures Figure 1.1 Predictions of final population size at four carrying capacities.................................24 Figure 1.2 Effect of data series length on predictions of final population size...........................25 Figure 2.1 IBM flowchart............................................................................................................57 Figure 2.2 Predicted pack composition.......................................................................................58 Figure 2.3 Predicted vagrant population size..............................................................................59 Figure 2.4 Predicted final population size...................................................................................60 Figure 2.5 Predicted hybridization rates......................................................................................61 Figure 2.6 Predicted probability of alpha pairs reproducing.......................................................62 Figure 2.7 Predicted final population size at various levels of management..............................63 Figure 2.8 Predicted final population size under various zone management scenarios..............64 Figure 2.9 Predicted probability of extinction at various levels of management........................65 Figure 2.10 Predicted hybridization rates at various levels of management...............................66 vii List of Abbreviations PVA Population viability analysis USFWS United States Fish and Wildlife Service ARNWR Alligator River National Wildlife Refuge IBM Individual-based model POM Pattern-oriented model ODD Overview, Design concepts, and Details (IBM model description) VSP Vagrant settlement probability MSP Mate selection probability VMP Vagrant movement probability MCP Management capture probability viii Chapter 1 : Effects of model structure and time series length on population viability predictions in the red wolf (Canis rufus) 1 Introduction Population viability analysis (PVA) can provide a powerful means of assessing viability and persistence potential within threatened or endangered populations (Brook et al. 2000, Beissinger and McCullough 2002, Morris and Doak 2002). However, uncertainties in parameter estimates and model structure associated with modeling stochastic environments can reduce the reliability of predictions made with PVA (Beissinger and Westphal 1998, Brook et al., 2000, Coulson et al. 2001). These uncertainties are often attributed to poor model specification, incomplete knowledge, and insufficient data (i.e. small sample sizes and short time series or data collection periods). As a results, the utility of PVA is often relegated to exercises in evaluating relative extinction risk within a competing management framework without any substantive contribution to the actual conservation objectives of the population. Yet, comparisons of relative risk may be insufficient as managers need to ensure that recovery objectives (i.e. minimum population size) are met in order to justify program costs, regardless of the relative merits of competing management actions (Brook et al. 2002). As more extensive demographic datasets are collected from small populations, more complete knowledge of threatened systems should facilitate greater accuracy in model structure and vital rate estimates, and in turn overcome many of the limitations of PVA suggested in the literature (Beissinger and Westphal 1998, Coulson et al. 2001, Doak et al. 2005). 1 However, conservation biologists often require accurate predictions of extinction dynamics prior to or during the early stages of data collection. In light of these needs, it is essential to evaluate the quality and accuracy of predictions made by initial PVA models developed during early periods in a conservation program's history. Brook et al. (2000) demonstrated the potential for predictive accuracy in PVA by estimating vital rates from partial datasets, while reserving the remaining data for an assessment of model performance. The results from the study, consisting of 21 demographic datasets from representative populations, suggest that short-term predictions of extinction risk and estimates of final population size are reasonably accurate across multiple systems, so long as the underlying parameters do not change substantially with time. In addition, population models often evolve to incorporate greater realism as more data become available. Further comparison between simple and complex models is needed in order to determine whether or not greater model complexity is even warranted within some systems. For example, Doak et al. (2005) demonstrated through simulation that deterministic models will generally outperform stochastic models with fewer than 5 years of data, provided one recognizes the overly optimistic tendencies of deterministic models. Although stochastic models provide more realism, they come with greater data requirements needed to accurately estimate the sources and extent of variation around vital rate means. However, incorporating stochastic processes only represents one way of changing model complexity. Changes in model structure, such as incorporating a social system or patterns of space use, often promote greater realism with substantial data requirement costs as well. Thus, an assessment of earlier, simpler models relative to more mature, data-driven models and observed trends within real populations may facilitate a better understanding of the efficacy of early PVA endeavors within specific systems. 2 Long-term monitoring programs associated with reintroductions often provide a wealth of demographic data that can be used to accurately parameterize population viability models (Gusset et al. 2009, Wakamiya and Roy 2009). Since the red wolf (Canis rufus) was first reintroduced into northeastern North Carolina in 1987 (USFWS 1989), the US Fish and Wildlife Service has closely monitored most of the individuals within the population through the use of VHF technology, traditional tags, and genetics. These efforts have contributed to an ever- growing demographic dataset encompassing a substantial period of time since the programs origins. In 1999, following almost 12 years of recovery efforts, participants in a population and habitat viability workshop generated the first substantial red wolf PVA (Kelly et al. 1999). Although not published within a peer-reviewed journal, the Kelly et al. modeling effort epitomizes what conservation managers often face when evaluating recovery objectives and management strategy, as well as population viability, during the early stages of management (i.e. short datasets, limited understanding of density dependence, etc.). The vital rate estimates used to parameterize their stochastic population model were correspondingly imprecise with broad confidence intervals that remained unchanged with changes in population density. Although the Kelly et al. model represented a preliminary assessment, there has since been no effort to evaluate the efficacy of their model in predicting population growth or viability in light of the limitations of the early red wolf dataset. With over 10 years of additional data since 1999, the red wolf reintroduction program provides a unique opportunity to build and parameterize a PVA with high quality data (i.e. longer data collection period, moderate sample size, almost complete demographic knowledge of the population, etc.). Here, we developed a simple population model that varies most substantively from the Kelly et al. model in the inclusion of a more biologically appropriate stage-structure and 3 density-dependent ceilings (Table 1). Our modeling effort was three-part. First, we re-evaluated the Kelly et al. (1999) PVA model performance relative to our more realistic and recently parameterized model, as well as to observed trends within the red wolf recovery area population, in an effort to assess differences in predictive outcomes based on modeling effort (i.e. personnel, model complexity, and data available at the time). Second, we evaluated the influence of study or data collection period length on viability predictions in the absence of an accurate, mechanistic model for density dependence within our more recently developed model. And third, we assessed the viability of the red wolf population in the recovery area with continued management. 2 Methods 2.1 Kelly et al. model description Kelly et al. (1999) used the population viability software package VORTEX (Miller and Lacy 2005) to predict population trends and evaluate the impact of wolf-coyote hybridization on red wolf population persistence. Although multiple models were constructed by Kelly et al., we chose to evaluate the only simplest model due to reservations about the structure of more complex models developed by the authors, which included hybridization between wolves and coyotes. In addition, an adaptive management plan had been implemented the following year to manage for wolf-coyote hybridization (Kelly 2000), potentially confounding an evaluation of model performance in the more complex models. The Kelly et al. model, henceforth the Workshop model, is described here to facilitate comparisons of model structure. The Workshop model represents a two-sex, stochastic population viability model incorporating two sources of environmental catastrophe: disease and hurricane. Vital rate estimates for the Workshop model were obtained from unpublished demographic data, professional knowledge, and published literature. Although Kelly et al. 4 estimated annual survival for both captive- and wild-born individuals, only survival estimates from wild-born wolves were used in their model. In addition, survival was stage-specific and was estimated for three appropriate age-classes: pups (0 ? 1 yrs), yearlings (1 ? 2 yrs), and adults (>2 yrs). Survival was subject to environmental stochasticity (unreported random distribution), but whether or not demographic stochasticity was included is unclear. The minimum and maximum breeding ages for red wolves were set at 2 and 8 years of age respectively. Only 50% (SD = 10% ; unreported random distribution) of females produced a litter within a given year. Similar to females, only 50% of males were arbitrarily assumed to be available for reproduction. Maximum litter size was set at 6 pups with an arbitrarily chosen set of probabilities for each possible litter size estimate (Table 2). Finally, the sex ratio of pups was assumed to be fixed at parity. Density dependence was incorporated in the Workshop model in the form of an absolute ceiling in population size (K = 150). When populations increased beyond this carrying capacity, additional mortalities were imposed equally across all age classes in order to return the population to K. In addition, two forms of catastrophe were incorporated in the Workshop model, represented by disease epidemics and hurricanes. Each event had an annual probability of occurrence of 1%. When either event occurred, survival and fecundity was reduced by 50% for that year only. 2.2 Current data and vital rate estimation The USFWS red wolf dataset is comprised of over 20 years of data (1987 ? 2007), encompassing the duration of red wolf recovery efforts in northeastern North Carolina. Since little is known about the demographics of red wolves prior to their reintroduction to the wild, all vital rate estimates used to parameterize the population models are potentially influenced by 5 management practices. Where possible, we controlled for management by for censoring known management influenced demographics (i.e. deaths due to management in survival analyses). Unless specified otherwise, all of the following models were parameterized using vital rate estimates from formal demographic analyses (Survival: Murray unpublished results, Reproduction: Steury unpublished results). The vital rates used to parameterize our model were estimated from 1987 ? 1993, 1987 ? 2000, and 1987 ? 2007 in order to evaluate the effect of time series length on predictions of viability (Table 2). Population size and trend estimates were used to compare predicted and expected population growth (USFWS 2007). These population size estimates did not represent rigorous statistical estimates of red wolf abundance, but rather minimum known red wolf counts during the fall of each year (USFWS 2007). Thus, we assume >80% of the population is known (Knowlton et al. unpublished results) and that there is no significant year-to-year bias in detection throughout the red wolf recovery area. 2.3 Current model description and process We generated a female-based, stage-structured population model, which incorporated demographic and environmental stochasticity in reproduction and survival with density- dependent social structure. Although Kelly et al. (1999) used VORTEX, we wanted the flexibility to explicitly define every major process within our model, as well as the ability to construct a simple social structure using a loosely individual-based approach. Therefore, we chose to use the statistical software programming language R ver. 2.9.2 (R Development Core Team 2004) for our model development. Since analyses of red wolf survival data suggest no difference exists in survival rates between the sexes (Murray unpublished results), and red wolf abundance is not sex-biased in the Red Wolf Recovery Area (Knowlton et al. unpublished results), a single-sex modeling approach 6 was chosen. Although a single-sex model deviates from the Workshop model (Kelly et al. 1999), that model assumes a 50/50 sex ratio with non-sex specific survival. Therefore, in essence, the Workshop model could have been modeled with only females as we do here. We constructed a 5- stage model composed of pups, yearlings, breeding adults (alphas), non-breeding adult pack members (betas), and vagrants (i.e. non-territorial floaters) in order to more accurately reflect social structure and pack composition within red wolves. A second model was built around red wolf age structure. However, results were very similar between age- and stage-structured models; hence, the more parsimonious stage-structured model is presented here. Annual survival was estimated by stage and birth location (i.e. captive-born, island-born, or wild-born; Table 2; Murray, unpublished results). As with the Workshop model, only wild- born estimates of survival were used in our model. The first process executed in the model schedule was survival in order to determine the number of surviving individuals used in later population/model processes. Stage-specific survival incorporated environmental and demographic stochasticity. For every time step, survival probabilities were drawn from beta distributions using stage-specific means and beta-adjusted standard deviations to impose environmental stochasticity (Morris and Doak 2002). Demographic stochasticity was integrated into the model by evaluating survival of each 'individual' in a Bernoulli trial framework using the environmentally-stochastic survival probability (Morris and Doak 2002). Our model allowed only surviving alpha females from the previous year to breed. Fecundity was modeled via three components: the probability of an alpha female reproducing, the number of pup recruits per litter in the fall, and the proportion of female pups. The annual probability of an alpha female reproducing was subject to environmental stochasticity and randomly determined at each time step using a beta distribution (Table 2). The annual probability 7 was in turn evaluated for all reproducing individuals using a Bernoulli trial framework to account for demographic stochasticity. The number of pup recruits was determined by drawing a random number of pups for each alpha female (demographic stochasticity; Table 2) from the most appropriate distribution parameterized using estimates determined from analyses of red wolf data (environmental stochasticity; USFWS unpublished results). In order to determine the number of pups that were female, the model once again used a beta distribution (Table 2) to specify the annual probability of being female (environmental stochasticity), which was in turn used to assess each pup in a Bernoulli trial (demographic stochasticity). The total number of logical positives represented the number of females recruits at each time step. Density dependence was incorporated in the model in order to provide a realistic limit on population growth. Demographic analyses have been unable to detect density dependence in survival or reproduction (Murray unpublished results, Steury unpublished results.). However, the extant red wolf population appears to be at carrying capacity with approximately 85-90 adults (Knowlton et al. unpublished results) in approximately 20 packs (estimated 115 ? 125 total individuals including pups, USFWS 2007; Analytically derived 138.7 individuals; Murray unpublished results). Since the mechanism for density dependence in red wolves has yet to be elucidated, density-dependence could not be explicitly modeled. With territorial species such as wolves, population modelers often impose limits (i.e. ceilings) on the number of territories within a population and attain realistic population dynamics (Gaona et al. 1998, Morris and Doak 2002). Therefore, two ceilings were applied to the social or stage structure of wolves within the model. The first was imposed directly by limiting the number of alpha or breeding females, thereby limiting the number of simulated packs. Four limits to the number of alpha females (i.e. 15, 20, 25, and 30) were implemented and the predicted population dynamics under each 8 imposed limit were compared with actual red wolf population trends in order to assess model performance (Fig. 1). A second ceiling explicitly modeled the number of betas within a pack, preventing unrealistically large pack sizes within the model, by randomly drawing a beta pack size number for each alpha female (i.e. pack) from an appropriate distribution (i.e. Poisson, zero- inflated Poisson, or zero-inflated negative binomial); the most appropriate distribution and parameters for the distribution were estimated from the red wolf data (Table 2). Incorporating pack composition and a more specific stage structure (pup ?12mo, yearling >12mo and ?24mo, and vagrant, beta, and alpha adults >24mo) represent the most significant deviations from the VORTEX model generated by Kelly et al. (1999). Pack size and composition was explicitly modeled in order to track stage assignments (i.e. social structure) and prevent unrealistically large populations. Pack composition followed a four step process. First, alpha individuals that died were replaced by non-alpha adults (i.e. yearlings, betas, and vagrants surviving from the previous year) within the density-dependent limits imposed by the model (i.e. 15, 20, 25, or 30 alphas). Second, if additional non-alpha adults remained, the number of alphas (i.e. packs) was used to generate the number of potential beta positions (i.e. # of beta adults per pack) available (see above). The number of betas per pack was simply halved to represent the number of available female beta positions (non-whole numbers were randomly rounded up or down using a Bernoulli probability = 0.5). Once the number of beta positions was determined, positions were filled by randomly selecting individuals from the pool of non-alpha adults. Third, all remaining non-alpha adults were assigned vagrant status. Fourth, all surviving pups transitioned into yearlings and were assumed to remain in their natal packs. The method used for model initialization depended on the conditions being evaluated. In general, however, 20 individuals were equally divided into 4 of the 5 stage classes (i.e. 5 9 individuals within pup, yearling, beta, and alpha stages). For comparison with the Workshop model, we initialized the simulations using 9 females divided equal into three stages: pup, beta, and alpha. A brief assessment of ergodicity was performed to assess the influence of initial conditions (i.e. initial population size and distribution) on model performance and system behavior. Once initialized, the model was advanced one time step (i.e. year) and the life history events discussed above scheduled. 2.4 Assessing model performance In order to capture the influence of stochasticity on population persistence, our model was iterated 1000 times over 25 and 50 year projections using annual time-steps. A quasi-extinction threshold of 2 total individuals from any stage was used to evaluate extinction dynamics and population risk within the simulated red wolf population. Our model was used to estimate viability statistics such as stochastic population growth (?s), mean final population size and variance, stage distributions, mean time to quasi-extinction, and the probability of quasi- extinction (i.e. proportion of iterations ending in quasi-extinction). Finally, the elasticities of ?s and mean final population size to proportional changes in survival, reproductive, and beta pack size parameters were evaluated using simulation (i.e. manual perturbation analysis) in a stochastic model framework and under density-independent and density-dependent conditions (Caswell 2001, Morris and Doak 2002). 3 Results 3.1 Model structure Before comparing performance in the two structurally different population models, we calibrated our model to the appropriate carrying capacity (K) or number of breeding females (i.e. packs). We ran our model allowing for 15, 20, 25, and 30 breeding females and compared model 10 output with observed population trends in the recovery area from 1990 ? 2006 (Fig. 1). A cap of 20 breeding females fit actual population trends best and reflects estimates for actual numbers of red wolf packs in northeastern North Carolina (USFWS 2007, Knowlton et al. unpublished results). Both models produced almost identical predictions of strong positive growth within the population (Workshop ?S = 1.215 ; Mahoney et al. ?S = 1.212) (Table 3), with our model predicting a more precise estimate of ?S in part due to greater precision in stochastic model parameters and vital rate estimates (Table 2). The ?S for both models were substantially lower than the reported observed population growth (? = 1.413 or rmax = 0.346; USFWS 2007). In truth, animals introduced into the population through management, especially during the early years of the recovery program, will artificially inflate any estimate of true population growth due to the physical presence and reproductive contributions of these individuals. Yet, the strong estimates of stochastic population growth for both models led to similar predictions of red wolf persistence within the recovery area over the next 50 years under continued management and in the absence of any major natural disaster (Prob. of quasi-extinction < 0.00 at Nquasi-extinction ? 2). As mentioned above, both models accurately reflect the observed stochastic population growth rate below carrying capacity and therefore follow observed density-independent trends quite well (i.e. 1990 - 1998). However, the differences in how density dependence was incorporated in the respective models led to the most significant deviation in model performance. The Workshop model imposed a limit of 150 total individuals (i.e. males and females of all stages). Not surprisingly, due to the arbitrary nature of this ceiling on total individuals, the Workshop model predicted a mean final population size at carrying capacity of 145 (SD = 13) individuals, 26 - 17% above the estimated 115 ? 125 total individuals within the actual red wolf 11 population. Our model predicted approximately 60 individual females and assuming a population in sexual parity, 120 total individuals within the population at carrying capacity (Table 3). 3.2 Length of data collection period In addition to evaluating the influence of model building effort and structure on viability predictions, we used our model to evaluate the predictive accuracy given vital rates estimated from three different data collection period lengths (i.e. 7, 14, and 21 years). Although our demographic analyses were unable to detect density-dependent effects on survival, the survival estimates across the three data period lengths suggested the presence of such an effect (drop in survival with for most stages with more data at or near carrying capacity; Table 2). In general, model predictions for mean final population size, variation in mean final population size, and population growth (?S) decreased with increasing time series length (Table 3). Population trajectories were accurate during the time periods over which vital rates were estimated (Fig. 2). The main difference in model performance between time-series lengths was highlighted in the predictions of projected final population size (i.e. population size predictions beyond the time period used to parameterize the model). The shorter time series appeared to bias predictions of population size relative to observed population growth over the associated time period (i.e. stronger or weaker growth lead to larger or smaller predicted population sizes, respectively; Fig. 2 a, b, and c). Not surprisingly, the model parameterized using the complete 21 years of data most accurately reflected short-term growth and observed trends within the red wolf population (Fig. 2d). All three time series used in the model predicted similar probabilities of quasi-extinction (Pquasi-extinction < 0.00 at Nquasi-extinction ? 2) (Table 3). Although study length made little difference in predictions of quasi-extinction in this red wolf population, it should be fairly apparent how larger 12 estimates of stochastic growth and final population size associated with the shorter time periods at low densities could lead to overly optimistic predictions of population persistence in other systems or populations (Fig. 2 A and B). 3.3 Perturbation analysis: the complete model The structure of our model lent itself best to manual perturbation in order to assess the relative influence of model parameters on population growth and final population size. For this analysis, we used our model parameterized with vital rate estimates from the complete data set (1987 - 2007). Table 4 outlines the elasticities of lambda (?S) and final population size (NFinal) to proportional changes in a single model parameter while holding all others to their predefined stochastic conditions. Following traditional approaches (Caswell 2001), we first evaluated the elasticities of lambda in the absence of density-dependence. Breeding adult survival (i.e. alpha females) appeared to have the greatest influence on population growth under density-independent conditions. In addition, pup and yearling survival, probability of an alpha female breeding, probability of a recruit being female, and the mean number of recruits all appear to be of equal importance, but secondary to adult survival. Although perhaps not as meaningful, we also evaluated the elasticity of ?S under density- dependent conditions (i.e. K = 20 breeding females) in order to identify potential changes in parameter importance relative to density-independent conditions. Survival of alpha females still appeared to have the greatest influence on population growth. However, the most notable differences between density-dependent and independent conditions are the relative reduction in the influence of yearling survival and increase in the influence of the probability of an alpha female breeding on ?S (Table 4), suggesting reproductive parameters may be more influential 13 than pup and yearling survival at higher densities. Since model predictions of final population size varied the most between the three time- series lengths, we evaluated the elasticity of NFinal to perturbations in model parameters (Table 4) under stochastic, density-dependent model conditions. Not surprisingly, the imposed limit on the number of breeding females (K) appears to have the greatest influence on estimates of final population size. However, perhaps more interestingly, the reproductive parameters for the probability of an alpha breeding, the probability of a recruit being female, and the mean number of recruits (in order of importance) all appear to substantially influence estimates of final population size under density-dependent conditions. 4 Discussion 4.1 Model Structure The Workshop model and our more recent model varied in complexity, programming language, and vital rate estimates used to parameterize the respective models. Although these represent several sources of variation, one of our objectives was to compare model performance and viability predictions given these pooled differences (i.e. compare predictions given different data, model developers, and software used). The most notable outcome was the almost identical predictions of stochastic population growth (?S) (Table 3). In addition, both models fit observed trends in the red wolf population over the time period during which vital rates were estimated. However, the Workshop model overestimated the mean population size relative to the observed population size at carrying capacity. This has almost entirely to do with the artificially-imposed ceiling on total individuals within the population with no biological feedback in the Workshop model. Incorporating a more biologically realistic and easier to estimate limit on the number of packs by restricting the number of breeding females, while modeling beta pack size explicitly, 14 more accurately reflects the observed population size at carrying capacity with 115 ? 125 total individuals (i.e. males and females of all stages). The results from this comparison demonstrate that remarkably similar conclusions can arise even provided substantial differences in model structure and logistical conception between two red wolf population models. Had there been greater differences in model predictions, an additional assessment would have been warranted in an effort to elucidate the cause(s) of any observed deviations (i.e. model structure, dataset, and/or modeling software). However, in light of our findings, we felt it was unnecessary, if not informative, to identify reasons for the similarities. If our goal was simply to evaluate persistence potential of the current red wolf population in northeastern North Carolina under continued management, than the simpler Workshop model would be sufficient and lead to similar predictions of viability. However, we believe our model will more accurately reflect system dynamics under more general scenarios (i.e. stochastic events or under conditions where such strong positive growth is not observed), especially if a more accurate sub-model for density-dependent survival can be incorporated. Until the latter has been accomplished, our model most realistically incorporates density dependence by limiting the number of reproductive individuals, and in turn packs within the population (an easier parameter to estimate than total individuals at carrying capacity). This simple method of imposing density dependence will allow for more accurate estimates of absolute population size following recovery area expansion or during future reintroduction efforts. 4.2 Length of data collection period Predictions of population growth (?S) and probability of quasi-extinction were in very close agreement across all time series model simulations, suggesting in the case of the red wolf, 15 time-series length may have little influence on predictions of population growth and persistence potential. However, our predictions of final population size for the first two time series were overly optimistic relative to observed population trends and likely due to our inability to mechanistically incorporate density-dependence within our model. As has been suggested in previous studies (Beissinger and Westphal 1998), accurately incorporating density dependence will be essential in developing more accurate estimates of final population size. Although we were unable to detect the influence of density on vital rates at the pack, regional, or population level (Murray unpublished results, Steury unpublished results), vital rate estimates from the four time series (three of different lengths) strongly suggest density- dependent survival in red wolves (Table 2). In fact, we would not expect to be able detect density dependence over the shortest time series (i.e. 1987 ? 1993, 1987 ? 2000, and 2001 - 2007) given the limited density-influenced survival data during these periods, representing one of the hazards of performing PVA with short or limited time series data. Although we were unable to determine the underlying mechanistic nature of population regulation in red wolves with the complete data set, and in turn a cause for the disparity between time-series estimates and preliminary, albeit formal, survival analyses, we propose several possible explanations: 1) an unmeasured area- or habitat-specific characteristic induces population regulation in red wolves, 2) weak density dependence in reproduction and survival that is undetectable in a single life history analysis alone, 3) management action occurs disproportionately more frequent at higher densities, and/or 4) poor model specification. However, until a more accurate, mechanistic sub-model for density- dependent survival is developed, our current model will be limited to evaluating population persistence within the current Red Wolf Recovery Area under imposed density-dependent limits. In scenarios where density dependence isn't modeled explicitly, we would caution that 16 modelers give careful consideration to time series length and where the population is on the density-dependent growth curve when estimating vital rates used to parameterize PVAs. In our case, the two shorter time series (i.e. 7 and 14 years of data) produced accurate population trends and estimates of population growth (?S) for the time period over which they were estimated, but predicted estimates of final population size overshot observed population carrying capacity, likely due to overly optimistic vital rate estimates associated with low density, exponential-like growth. Intuitively, vital rates estimated from shorter time series should strongly reflect population trends over the time period during which the data was collected. In order to demonstrate this interaction, we re-estimated our vital rates over a fourth time series from the years 2001 through 2007 (assuming little knowledge of or data from previous years) during which red wolves have experienced a decline. As one can see in Fig. 2C, our model using vital rates estimated from the last seven years predicted substantially lower population growth and final population size. In addition, although our model parameterized with vital rates estimated from 21 years of data accurately predicted population trends and reflected current estimates of carrying capacity in red wolves, the complete model slightly underestimated realized population growth at population sizes below carrying capacity (Fig 2C). Thus, stochastic vital rates estimated from populations currently at carrying capacity may in fact negatively bias estimates of population growth and in turn increase predictions of extinction risk following disturbance. The results from the perturbation analysis of our complete model suggest population growth (?S) is most elastic to changes in adult survival under density-independent conditions. However, the red wolf population in northeastern North Carolina is currently at carrying capacity, suggesting other vital rates may be important in maintaining population levels. Although not commonly done, we evaluated the elasticity of ?S and NFinal to perturbations in 17 model parameters under density-dependent conditions (K = 20). The results from the latter two analyses were in agreement, suggesting sustainment of strong reproduction will encourage strong growth on average and larger population sizes. Interestingly, one would expect vagrant survival to have a larger effect on estimates of population growth and final population size given vagrants are not subjected to imposed density limits. We suspect the discrepancy is due to the limited production of, and to some extent low survival for, vagrants within the model. Reproductive parameters likely become more important at carrying capacity due to greater production of surplus animals, and in turn vagrants, inflating estimates of final population size. Therefore, maintaining strong reproduction at carrying capacity may be not only important for maintaining genetic diversity and reducing the detrimental genetic effects of small populations by encouraging larger population sizes, but also in reducing the probability of hybridizing with non- wolves by increasing regional densities and the availability of mates. In conclusion, differences in model structure and length of time series used in estimating model parameters would ultimately have little influence on our predictions of red wolf persistence in northeastern North Carolina. In general, there was strong agreement among all models for predictions of stochastic population growth (?S), perhaps the most commonly used and supported measure of population viability. Where models deviated the most were in predictions of mean final population. In a world where absolute benchmarks towards recovery are often needed (i.e. minimum population sizes, etc.), predictions of final population size are a desirable outcome of PVA. The differences in predictions of mean final population size are attributed to our inability to explicitly model density dependence within the system. First, how we impose density dependence (i.e. model structure) can lead to very different predictions of absolute population size. The Workshop model and our more recent model using the complete 18 data set produced similar estimates of final population size, likely due to strong a priori knowledge of carrying capacities within the red wolf recovery area. With incomplete or incorrect knowledge of the system, however, the results could have been drastically different, stressing the importance of having a strong understanding of the carrying capacities or biologically-limiting 'controls' within the population(s) before efforts to impose density dependence are made (i.e. absolute ceilings on total individuals or specific stages). Second, without explicitly modeling density dependence, vital rate estimates can be strongly influenced by time series length in a density-dependent population. As in our case, shorter time series can lead to overly optimistic or conservative vital rate estimates when population densities are changing, which in turn can lead to biased estimates of mean final population size. The results from this modeling endeavor once again stress the need to explicitly model density dependence in viability models for appropriate populations. If one is unable to do so, as is often the case, than the potential biases must be acknowledged and/or addressed. 19 Table 1: An outline of key similarities and differences between two red wolf population viability models. P, Y, A, B, NB, and V represent pup, yearling, breeder, non-breeder, and vagrant stages respectively. NTotal - total population size. NB - number of breeding females (alphas). NBetaPack - number of non-breeding pack members (betas). 20 Table 1 : Structure Schedule Post-breeding Post-breeding Sex M/F F Carrying Capacity Stochasticity Environmental Environmental/Demographic Data Years 1987 ? 1998 1987 ? 2007 Software VORTEX R Kelly et al. Workshop Model Mahoney et al. Model Age/Stage (P,Y,A) Stage (P,Y,B,NB,V) NTotal = 150 NB = c(15,20,25,30)N BetaPack ~ ZI Pois or ZI NegBinom Table 2: Vital rate estimates used to parameterize Kelly et al.'s baseline model (workshop model) and our stage-structured model with associated estimates of error (SE). N represents the sample size used in individual analyses. aWhereas the workshop model estimated the proportion of all females that breed, our estimate is the proportion of alpha females that breed. bParameter estimates for the number of recruits in the Fall. N/E ? Not estimated. Pois ? Poisson distribution parameters. ZI Pois ? zero-inflated Poisson distribution parameters. ZI Neg. Binom ? zero- inflated negative binomial distribution parameters. P0 ? probability of zero. 21 Table 2 : Survival (Wild-born) 1987 - 1997 1987 - 1993 1987 - 2000 1987 - 2007 Adult 0.87 (0.10) N/E N/E N/E Adult, dominant N/E 0.952 (0.042) 0.947 (0.032) 0.901 (0.032) Adult, beta N/E 0.553 (0.076) 0.780 (0.042) 0.649 (0.020) Adult, vagrant N/E Yearling 0.83 (0.10) 0.926 (0.053) 0.840 (0.034) 0.806 (0.026) Cub 0.78 (0.10) 0.881 (0.056) 0.897 (0.028) 0.816 (0.041) N = 87 N = 274 N = 435 Reproduction Prob. of female Breeding 0.500 (0.100) Prob. of a female cub 0.500 (0.100) Number of pups Pack Size Beta pack size N/E Kelly et al. Workshop Model Mahoney et al. Model 0.006 (3.435e-06) 0.222 (5.440e-04) 0.164 (1.266e-04) 0.822 (0.290)a N = 20 0.836 (0.186)a N = 106 0.843 (0.026)a N = 197 0.565 (0.267) N = 42 0.472 (0.188) N = 212 0.457 (0.015) N = 501 1 (0%), 2 (20%), 3 (20%), 4 (40%), 5 (50%), 6 (10%) Poisb ? = 2.563 N = 16 ZI Poisb ? = 2.464 P0 = 0.074 N=88 ZI Neg. Binomb 2.375 (k=8.579) P0= 0.248 N = 185 Pois ? = 0.773 N = 44 ZI Pois ? = 1.201 P0= 0.309 N = 219 ZI Pois ? = 1.197 P0= 0.425 N = 488 Table 3: Model predictions by model type and time series length. Mean NFinal is the mean final population size and either includes males and femalesb or just femalesf. ?S is the mean stochastic population growth rate. All Mahoney et al. model means estimated from 1000 iterations. 22 Table 3 : Data Years 1987 ? 1997 1987 ? 1993 1987 ? 2000 1987 - 2007 1.22 1.23 (0.037) 1.226 (0.024) 1.212 (0.028) <0.00 <0.00 <0.00 <0.00 Kelly et al. Workshop Model Mahoney et al. Model Mean NFinal (SD) 145 (13)b 87.166 (24.375)f 77.182 (15.146)f 59.892 (11.803)f ?S (SD) Prob. Of Quasi- extinction (N ? 2) Table 4: Results from the manual perturbation analysis of the complete Mahoney et al. model (data from 1987 - 2007). Values represent elasticities for ?S and NFinal to proportional changes in single model parameters. 23 Table 4 : Survival Parameters Pup Survival 0.24 0.04 27.51 Yearling Survival 0.23 0.01 18.75 Vagrant Survival 0 0 -0.63 Beta Survival 0.01 0.01 7.07 Alpha Survival 0.63 0.1 29.91 Reproductive Parameters Prob. of alpha breeding 0.26 0.07 44.17 Prob. of female recruit 0.24 0.04 41.22 0.24 0.04 39.93 Size (ZI Neg. Binomial) -0.01 0.01 -3.43 -0.07 -0.04 -9.07 Pack Size Parameters 0 0 8.05 Prob. Of Zero (ZI Poisson) -0.01 0.01 3.61 K N/E N/E 59.79 Density- independent Lambda Density- dependent Lambda Density- dependent NFinal Mean # of recruits (ZI Neg. Binomial) Prob. of zero (ZI Neg. Binomial) ? (ZI Poisson) Figure 1: Mahoney et al. model simulations (1000 iterations) at K = c(15, 20, 25, 30) breeding females using vital rates estimated from the complete data set (1987 ? 2007). The solid lines represent predicted mean population size at each year in the simulation with dashed lines at one standard deviation. The solid black line shows the estimated number of females within the actual red wolf population (50% total estimates; USFWS 2007). 24 Figure 2: Mahoney et al. model simulations (1000 iterations) at K = 20 breeding females using vital rates estimated from three time series lengths (four time periods). The solid lines represent predicted mean population size at each year in the simulation with dashed lines at one standard deviation. The dotted line shows the estimated number of females within the actual red wolf population (50% total estimates; USFWS 2007). 25 Chapter 2: The role of behavior and management in predicting red wolf (Canis rufus) persistence potential: an individual-based modeling approach 1 Introduction Classical population viability (PVA) and life history analyses often rely on system-level processes, such as average survival and reproduction, to define overall population dynamics, essentially fitting models to observed trends within populations (Beissinger and McCullough 2002). Well designed models should replicate observed temporal data, but will not necessarily provide reasonable predictive accuracy (Beissinger and Westphal 1998). The inaccuracy of population viability models is in part due to the often very specific conditions in which data were collected, potentially leading to poor predictive power and reduced generality across systems (Burgman and Possingham 2000, Possingham et al. 2002, Wiegand et al. 2003). In addition, classical models that focus on system-level processes may ignore patterns within individuals that are critical to population persistence (Vucetich and Creel 1999, Grimm and Wissel 2004, Gussett et al. 2009). Consequently, such models may have limited applicability to real systems (Caughley 1994, Wiegland et al. 2003). Individual-based models (IBM), or pattern-oriented models (POM), overcome some of the limitations of classical PVA (Grimm and Railsback 2005). Modeling ecological systems by reproducing patterns within individuals in an IBM framework has several notable advantages. First, the framework of an IBM is built from the ?bottom-up? allowing one to incorporate and observe how patterns in behavior influence emergent, population-level dynamics with limited 26 imposed system processes (Grimm and Railsback 2005). Second, by accurately modeling patterns in individual behavior, one can hypothetically eliminate or reduce system-specific elements from the model, thereby increasing the generality, applicability, and realism of the model. Third, IBMs appear to be robust to parameter uncertainty, potentially overcoming model limitations associated with deviations in demographics between systems (Grimm and Railsback 2005). In light of these traits, IBMs can be powerful tools for evaluating competing management strategies and population performance within endangered species where more general models are needed (Wiegand et al. 1998, Lacy 2000, Gusset et al. 2009, Swanack et al. 2009), such as in site-specific reintroductions of rare species (Kramer-Schadt et al. 2004, Gusset et al. 2009). The experimental population of red wolves (Canis rufus) in northeastern North Carolina is an excellent example of a population with numerous complexities at the level of the individual that are integral to an evaluation of population performance and would be best studied using an IBM. The red wolf (Canis rufus) historically occupied much of the deciduous landscape of the eastern United States (Nowak 2002). Due to intensive predator control efforts by humans following colonial settlement of North America, red wolves were reduced to small remnant populations along the gulf coast of Texas and Louisiana by the late 1960s. Within six years of having been listed as endangered in 1967, US Fish and Wildlife Service biologists concluded that red wolf persistence in the wild was unlikely due to hybridization with coyotes (C. latrans) and initiated an intensive capture effort of all suspected red wolves throughout the southeastern US (USFWS 1990). Following the successful establishment of a captive breeding program, the US Fish and Wildlife Service reintroduced the first experimental population of red wolves in Alligator River National Wildlife Refuge, North Carolina in 1987 (USFWS 1989). The red wolf population has since expanded throughout 5 counties within the Albemarle Peninsula of 27 northeastern North Carolina (~6,000 km2) to an estimated 85-90 adult red wolves in approximately 20 packs as of 2009 (U.S. Fish and Wildlife Service 2007, Knowlton et al. unpublished results). However, with increases in coyote abundance through the mid-1990's, the historical problems associated with coyotes returned (Phillips et al. 2003, Adams et al. 2007), and by 1999 coyote hybridization and introgression was recognized as the most significant threat to red wolf recovery efforts (Kelly et al. 1999). In response, the USFWS developed an adaptive management plan to mitigate the impacts of hybridization and introgression on red wolves through zone-specific sterilization or euthanasia of coyotes within the recovery area (Kelly 2000; Currently: Zone 1 euthanasia ; Zones 2 & 3 sterilization). The adaptive management plan remains active and has proven effective at reducing the impacts of hybridization and spatial competition with coyotes (Fredrickson and Hedrick 2006). Recent population-level assessments of red wolf viability demonstrate strong growth below carrying capacity within the experimental population of red wolves (Kelly et al. 1999, Mahoney et al. in prep). However, due to limitations of classical viability approaches, previous models were unable to accurately account for complex socio-behavioral patterns that are likely important for population persistence in the red wolf. Hence, a more general model that incorporates the behavioral dynamics of this two-species system is needed to better understand how red wolves will persist with coyotes on the landscape, as well as to meet the additional recovery objectives outlined in the Red Wolf Species Survival Plan (USFWS 1990). Recently, individual-based models have been developed to model territoriality and social structure within canids (Vecetich and Creel 1999, Pitt et al. 2003, Gusset et al. 2009), as well as to explore alternative predator control strategies for coyotes (Conner et al. 2008). In addition, individually- based methods have been used to explore competitive space-use in the red wolf-coyote system 28 (Roth et al. 2008). Models such as these can be expanded to incorporate system-specific behavioral elements (i.e. mate choice and space use) within a competing canid species framework (i.e. red wolves and coyotes), and implemented as an evaluation of population viability for one or both species. Behavioral processes such as space use and dispersal, social behavior, and mate choice may be difficult, if not impossible, to incorporate accurately using classical population models, but likely play a significant role in population persistence (Kelly et al. 1999, Grimm et al. 2003, Gusset et al. 2009). Furthermore, various management actions associated with in the red wolf recovery program have effects on both coyote and red wolf individuals that can only be fully understood via an examination of how such programs influence behavior, specifically space use and social interactions with other canids. Therefore, we developed an individual-based model that explores the impact of complex socio-behavioral patterns, as well as management effort, on red wolf viability in the presence of coyotes. Although some of these behaviors may be unique to this system, exploring individual behavior within the framework of IBM can provide further insight into how social interactions impact general population-level processes. 2 The Model The model description provided below follows the ODD protocol (i.e. Overview, Design concepts and Details) recommended by Grimm et al. (2006, 2010) for describing individual- based models. The main IBM was built using the pre-existing agent-based modeling platforms NetLogo ver. 4.1.1 (Wilensky 1999) and NetLogo-R-extension (Thiele and Grimm 2010) in order to expediate model building, as well as to provide a more efficient means of interpreting model outputs. All of the supporting analyses and sub-model building were completed using the statistical software R ver. 2.11.1 (R Development Core Team 2004). 29 2.1 Purpose of model In building the IBM, our goal was to evaluate red wolf population viability in the context of management by incorporating individual space use and socio-behavioral patterns from the coyote and red wolf system. We were interested in behaviors associated with pack structure/formation and mate replacement, as well as how each in turn was influenced by management, and their respective roles in red wolf population dynamics. More specifically, we used the IBM calibrated to recovery area dynamics in order to evaluate red wolf population growth and extinction risk under various management regimes. We further evaluated how probability of non-wolf capture and seasonal capture effort influenced population performance in red wolves. We also explored alternative zone management scenarios by examining various spatially explicit coyote control efforts. Finally, we ran the model under three mate replacement regimes (i.e. random, semi-assortative, and full assortative) in order to assess the spectrum of variability in community dynamics associated with hybrid pairing and its interactions with various management protocols. 2.2 Entities, state variables, and spatial scales Three hierarchical entities were established and monitored within the IBM: individual canids, packs, and populations (all three were species-specific and spatially-defined). The primary fundamental unit was the individual canid, which was characterized by the following state variables: individual ID, species (coyote, red wolf, or hybrid), sex, age (in months), stage (alpha ? dominant pack member, beta ? non-dominant pack member, yearlings ? between 12 and 24 months, pups ? <12 months, and vagrants; all mutually exclusive), number of prior dispersals (influences dispersal probability ? see below), time since last dispersal (also influences dispersal probability), breeding condition (sterile or intact), current pack membership, and whether or not 30 the canid was in its natal pack (for pack relatedness and beta proportion estimates). Only red wolves ?12 months old could attain alpha status (yearlings attaining such status were no longer considered yearlings). Coyotes were defined as pups at ages ?6 months and transitioned to adult stages at >6 months following Conner et al. (2008). Thus, coyotes could attain alpha status as adults much sooner than red wolves. Since the genetic composition of individuals was not explicitly defined within the model, all hybrids were treated as 'coyotes' by assigning coyote- specific state variables. The 'pack' represented the fundamental reproductive unit ? each pack could have only one alpha male and one alpha female and only alphas were allowed to reproduce. Each pack was treated as a single entity and assigned the following state variables: pack ID, pack type (wolf ? both alphas were wolves, coyote ? both alphas were coyote, or hybrid), location, list of pack members, presence and ID of an alpha female, presence and ID of an alpha male, ID of the alpha female in January, and breeding status (was the alpha female impregnated and capable of having a litter of pups). Space was represented by a 5x6 grid of square cells, with each cell assumed to represent an average pack territory. An additional 30 territories were incorporated along one side to serve as a dispersal buffer zone for wolves and population source for coyotes (beyond model initialization, no new coyote individuals were spawned other than those produced during reproduction). This particular spatial arrangement was chosen in order to loosely replicate conditions within the Red Wolf Recovery Area (i.e. ~20 - 30 packs with immigration and emigration through a single western boundary; U.S. Fish and Wildlife Service 2007, Knowlton et al. unpublished results). Since the red wolf recovery area is effectively a peninsula, the three remaining sides of the assessment area were hard boundaries with no movement and/or dispersal 31 beyond these edges. Unless otherwise stated, model performance was assessed using the 'eastern' grid of 30 packs only. 2.3 Process overview and scheduling We used discrete intervals of one month in our IBM since most of the socio-behavioral characteristics modeled occur over shorter time scales than a year. The model year commenced in January. During one monthly time step, sub-models were implemented to evaluate the following processes (in order): individual survival, mate replacement (Alpha), modeler-imposed management action, reproduction, dispersal, vagrant movement, and pack establishment (Fig. 1). Packs were selected at random for processing, regardless of species makeup and without replacement until all packs had been assessed. Every individual was modeled from birth to death. 2.4 Design concepts 2.4.1 Emergent properties: The spatial and temporal dynamics within and between the two populations of canids emerged from explicitly-defined red wolf and coyote behavioral rules. In addition, behavior drove the emergence of species-specific population characteristics, such as pack composition, alpha pair composition, proportion of alpha pairs breeding, hybridization, total population carrying capacity, pack establishment rates, age and stage structure, and population persistence potential. 2.4.2 Sensing (individual knowledge): Individuals within a pack were assumed to have complete knowledge of neighboring territories, such as quality, species occupancy, pack size, and the social status of individuals (i.e. reproductive pair intact or disrupted). Such information was used by individuals in making 'decisions' about dispersal, movement, pack establishment, and pack settlement. This assumption is likely realistic due to the chemical and physical cues used by wolves to demarcate boundaries, as well as from the exploratory movements often exhibited 32 prior to dispersal (Mech and Boitani 2003). 2.4.3 Individual interactions: No interactions between individual canids were defined explicitly in our model. However, numerous interactions could be accounted for implicitly. For example, displacement of an alpha, as well as the eviction of a beta, were subsumed in dispersal probabilities. Aggressive interactions between wolves were implicitly accounted for in mortality, dispersal, and vagrant movement probabilities. Mate replacement and competition for dominant status was subsumed into breeding pair formation. Alpha wolves implicitly suppress reproduction within subordinate pack members (i.e. non-alpha canids do not reproduce in the model). Finally, copulations were assumed to occur whenever an alpha female and alpha male were present during the month of January. 2.4.4 Stochasticity: The model contained no environmental stochasticity. All behavioral rules and demographic parameters, with the exception of the number of recruits, were based on probabilities (Table 1). Thus, demographic stochasticity for all probabilities was generated by using a random uniform distribution in a Bernoulli trial framework. Demographic stochasticity in the number of wolf recruits (pups) per breeding female was generated using a zero-inflated negative binomial distribution (see below; Table 1). 2.4.5 Observation: For model performance assessments, population structure and individual demographic traits were monitored and data collected during each time step from the core assessment area, with specific attention towards mate replacement, pack composition, and space use in the two species. 2.4.6 Initialization: The model simulations were initialized with 30 red wolf packs, or 30 alpha pairs, occupying the extent of the core assessment grid (i.e. 5x6 grid representing the recovery area). The buffer zone, cells occupying one side of the assessment grid, was initialized with 30 33 coyote packs. Additional pack members were randomly assigned using a zero-inflated Poisson distribution fitted to data on adult red wolf beta pack size (? = 1.197 , PZero = 0.425 ; Table 1). The sex of each additional individual was assigned at random assuming a population at sexual parity. Simulations were initialized during the month of December (January as the first time step) with randomly assigned ages reflective of an April birth (uniform distribution discretized by rounding to appropriate ages in months; 8mo, 20mo, 32mo, etc.). 2.4.7 Environmental input: Environmental variables were not included in the IBM. 2.5 Sub-models Unless specified otherwise, all sub-models were parameterized from analysis of data from the reintroduced red wolf populations. The USFWS red wolf dataset is comprised of over 20 years of data (1987 ? Present), encompassing the duration of red wolf recovery efforts in Alligator River National Wildlife Refuge, North Carolina. Since little is known about the demographics of red wolves prior to the involvement of management officials, all vital rate estimates derived from analyses using data from the experimental population are potentially influenced by management practices. 2.5.1 Mortality: A semi-parametric Cox proportional hazard analysis of red wolf data identified a significant effect of stage on mortality (Mahoney unpublished results). Therefore, we estimated monthly mortality rates from the data as a function of stage using Kaplan-Meier statistics, with the stage-specific estimates shown in Table 1 (Mahoney unpublished results). Mortality estimates for coyotes followed those used in previously published models (Pitt et al. 2003, Conner et al. 2008). Coyote mortality in adults (including vagrants) is a curvilinear function of age, while mortality in pups is fixed at 0.1 (Table 1). Mortality in vagrants is also a function of the number of vagrants and the total number of packs in the population (both wolves and coyotes; Table 1). 34 2.5.2 Alpha & mate replacement: The general process described here was the framework for all three mate replacement procedures used in the model. The difference between the three procedures was in how individuals were pooled at respective levels of the mate replacement process. The random procedure pooled individuals irrespective of species, the priority procedure gave red wolves priority with other red wolves only within a given mate pool, and the full assortative procedure only pooled individuals of the same species. Initially, replacement of open alpha positions gave priority to beta pack members of the same sex as dictated by the available position. Beta pack members (>24 mo for red wolves, >6 mo for coyotes) were ranked by age and individually given an opportunity to occupy the available position with a preset mate selection probability (MSP ; assessed using a random uniform distribution in a Bernoulli trial). If there were no betas or the betas were unsuccessful at changing social status, vagrants associated with the pack territory were ranked according to age and given an opportunity to fill the alpha position with the same probability, MSP. If there were no vagrants within the territory or they failed to occupy the alpha position, betas from neighboring territories were pooled, ranked, and 'tested' using the same method. If there were no betas available within the neighborhood (i.e. surrounding 8 territories) or they failed to change social status, vagrants associated with the neighborhood were given an opportunity to occupy the alpha position. Finally, if there were no other available canids or they failed to change status, yearling red wolves within the pack had an opportunity to fill the available alpha position. Although sex-biased hybridization has been suggested in a similar system (Canis latrans and Canis lycaon; Kays et al. 2010, Rutledge et al. 2010a), the model assumed equal probability of hybridization between the sexes as suggested in preliminary analyses for red wolves and coyotes in North Carolina (Mahoney et al. unpublished results). Since the mate selection probability (MSP) could not accurately be estimated from data 35 on red wolves or coyotes, we explored the effect of this nuisance parameter on model results (see below). 2.5.3 Management: The spatial core of the simulation area (eastern 30 cells) was divided into 3 management zones (Kelly 2000, U.S. Fish and Wildlife Service 2007). All non-wolf canids in the spatial core of the model (Zones 1 - 3) were assessed individually using a random uniform distribution with some probability of capture during every monthly time step within a management season. Those individuals that failed a Bernoulli trial were 'successfully' captured and killed in zone 1 or sterilized in zones 2 and 3 as currently implemented in the recovery area. 2.5.4 Reproduction: Litter production was only possible in packs with an intact breeding pair (i.e. 1 alpha male and 1 alpha female) during both the months of January and April. In addition, the alpha female in April must have been the alpha female in January in order for copulation to have occurred. For wolves, the number of recruits per pack rather than litter size was modeled as a response variable due to the higher quality data of the former. Reproductive analyses did not support a relationship between pack size and the number of fall recruits within red wolves (Steury et al. unpublished results). Therefore, we used a zero-inflated negative binomial distribution (the best fitting distribution; ?=2.375, k=8.579, and Pzero=0.248) to randomly determine the number of wolf recruits for each successfully reproducing pack (Table 1). Wolf pup recruits were generated in April, but were not susceptible to mortality or dispersal until the fall (in order to maintain consistency with fall recruit estimates). The sex of individual pups was determined using a random uniform distribution (i.e. Bernoulli trials with Pfemale = 0.50). Reproductive estimates for coyotes followed those used in previous models with the number of pups defined as a function of pack size (Pitt et al. 2003, Conner et al. 2008). Unlike wolf pup 'recruits', coyote pups were permitted to die, as well as disperse, prior to the fall. Reproduction 36 within 'hybrid' packs, specifically reproductive pairs composed of one red wolf and one coyote, were defined by the female's species. 2.5.5 Dispersal: Dispersal probabilities were incorporated in the IBM to reflect dispersal initiation of red wolves within the recovery area, as well as to enforce realistic values for mean red wolf pack size. The dispersal function was estimated with red wolf data from the recovery area using mixed logistic regression while allowing for multiple dispersals by a single individual (i.e. multivariate failure-time analyses; Mahoney unpublished results). Parameter estimates were averaged across all possible combinations of fixed effects, but always with a random effect of individual. Red wolf dispersal probabilities were estimated as a function of age, stage, pack size, number of previous dispersals, time since last dispersal, and season (Table 1). As with previous sub-models, coyote dispersal probabilities will follow the sub-model used by Conner et al. (2008) and were a function of pack size. 2.5.6 Vagrant movement: In order for alpha replacement by and management action on vagrants to be realistic, each vagrant was assigned a focal territory. During the vagrant movement action loop, vagrants were chosen at random and given an opportunity to stay or make a single move to a neighboring territory (i.e. the surrounding 8 territories) based on a system of priority rules. First, if a vagrant was within a pack-less territory with a canid of the opposite sex, the vagrant would remain in current territory. Second, the vagrant would move to any neighboring territory with an available alpha position of the same sex. Third, the vagrant would move to any neighboring territory unoccupied by a pack, but with a canid of the opposite sex. Fourth, the vagrant would move to any unoccupied neighboring territory. And fifth, the vagrant would randomly move to any neighboring territory with a preset probability regardless of occupancy (Table 1). Because this vagrant movement probability (VMP) could not be estimated from data, 37 we explored the effect of this nuisance parameter on model results (see below). 2.5.7 Pack establishment and vagrant settlement: If two vagrants of the opposite sex were associated with an empty focal territory (i.e. no pack present), a new pack was established regardless of species (except models with full assortative mate replacement). Pair formation during pack establishment followed the same priority rules as specified by the pair formation procedure during alpha replacement (2.5.2). In addition, potential red wolf pairs were given priority over coyotes for space due to the larger size and dominant status of the former species (Fredrickson and Hedrick 2006, Rutledge et al. 2010a). If more than two vagrants of the same sex and species were present, alpha status was randomly assigned to one while the remaining individuals kept their vagrant status. Finally, vagrants within territories occupied by packs were given a preset probability of joining the resident pack (Vagrant Settlement Probability, VSP). Once again, since this vagrant settlement probability (VMP) could not be estimated from data on red wolves or coyotes, we explored the effect of this nuisance parameter on model results (see below). 3 Simulations & Results 3.1 Model robustness to uncertainty in nuisance parameters and model structure In an effort to assess model robustness to unknown or inestimable parameters, BehaviorSpace in NetLogo was used to evaluate the sensitivity of several emergent properties to broad changes in three nuisance parameters (i.e. mate selection probability (MSP), vagrant movement probability (VMP), and vagrant settle probability(VSP)). Specifically, we evaluated the effects of variation in the nuisance parameters on mean red wolf pack size, mean number of betas per pack, mean number of yearlings per pack, mean number of red wolf vagrants, proportion of non-natal betas, total red wolf population size, mean number of new hybrid packs, 38 and the probability of a red wolf alpha females denning. Model simulations with random, red wolf priority, and full assortative mate replacement procedures were each evaluated over 200 iterations of 50 year projections for each parameter set with five levels of monthly coyote/hybrid capture probability (0.0, 0.25, 0.50, 0.75, 1.0) during the spring and fall seasons. Further, since model complexity and the quantity of model output limited our effort to a coarse assessment of parameter space, we regressed the three nuisance parameters on output from year 50 for each of the emergent properties using generalized linear models with appropriate error distributions. The regression models were then used to predict emergent property means at a finer scale for each mate replacement procedure, as well as to determine the proportion of model predictions within one standard deviation of observed means (an estimate of how robust the IBM was to inestimable parameters and model structure; Grimm and Railsback 2005, Kramer-Schadt et al. 2007). We found that variation in vagrant movement probability did not have a biologically significant effect on model performance for any emergent property regardless of mate replacement procedure (i.e. random, red wolf priority, and full assortative mating). Therefore, hereafter we limit our discussion to the impact of VSP and MSP on model performance. Additionally, we found that VSP and MSP showed consistent trends across the three mate replacement procedures for each emergent property examined. Thus, only results for the red wolf priority mate replacement procedure are depicted with a brief discussion of key differences between mate replacement regimes provided below. Mean red wolf pack size increased non-linearly with increases in both VSP and MSP, with larger effects at high and low probabilities respectively (Fig. 2a). The mean number of betas within a pack experienced similar trends with increases in VSP and MSP, although with a 39 substantially weaker effect of MSP (Fig. 2c). The number of yearlings within a pack demonstrated intermediate patterns with a larger effect due to MSP than was observed in beta pack size and a positive asymptotic trend in VSP (Fig. 2b). The number of pups within a pack was not evaluated since the production of pups was imposed by the model. Model predictions for pack composition at all three levels were very similar between priority and full assortative mate replacement procedures, with slightly smaller numbers on average under random mate replacement. The probability of a pack containing betas that were born to that pack (probability of natal betas) exhibited somewhat more complex behavior (Fig. 2d). Vagrant settlement probability had a negative, non-linear effect on the probability of natal betas with the largest effect at low to intermediate VSP. As with beta pack size, MSP had a small positive asymptotic effect on the probability of natal betas. Model predictions for probability of natal betas were almost identical between mate replacement procedures except at VSPs near one where an increasingly stronger positive asymptotic effect of MSP boosted natal beta proportions from random to full assortative mate replacement. Vagrant population size decreased non-linearly with increases in VSP (Fig. 3). In general, the number of vagrants increased asymptotically with an increase in MSP, but overall the effect of MSP decreased to almost non-existent as VSP approached 1.0. The total red wolf population size at 50 years was more complex with heteroschedastic data (Fig. 4 : the small hump at zero wolves represents populations that were extinct or going extinct at low management capture probability). In general, variability in final population size increased with increases in VSP. In addition, there were small non-linear increases in final population size with increases in VSP and MSP. Most notable, however, was the apparent interaction between VSP 40 and MSP, whereby increases in VSP disproportionately impacted population size at low MSP (i.e. 0.10 ? 0.25), resulting in more negative shifts in total wolves relative to moderate and high MSP (Fig. 4). Random and priority mate replacement procedures predicted almost identical final populations sizes at all levels of MSP and VSP. However, results for random and priority mate replacement procedures were lower than that under full assortative mate replacement (?20-30 individuals). Unlike the other emergent properties, there were very large differences in model predictions for hybridization rates (i.e the mean number of new hybrid packs per year) between random and red wolf priority mate replacement procedures (full assortative did not allow for hybrid pair formation). Although the overall trends between MSP, VSP, and hybridization were the same, with declines in hybridization as MSP increased and VSP decreased, the effect of MSP was much stronger relative to VSP under priority mate replacement and vice versa under priority mate replacement (Fig. 5). Finally, the mean probability of red wolf alpha pairs or packs reproducing was strongly influenced by both VSP and MSP (Fig. 6). An increase in VSP produced slightly non-linear declines in the probability of denning. However, increases in MSP asymptotically increased the probability of alphas denning with very strong rates of increase at low to moderate MSP. Here again, model predictions were very similar between all three mate replacement procedures. Although initially somewhat surprising, the random and priority mate replacement procedures essentially had the same 'total' mate pool with individuals simply re-arranged based on priority rules. With full assortative mate replacement, the red wolf population was able to maintain a sufficiently larger size to make up the difference coyote individuals would normally make in the mate pool and allow for comparable probabilities of reproducing. 41 Overall, the IBM was robust to uncertainty in model structure associated with mate replacement procedures, with only hybridization rates resulting in substantially different predicted effects of MSP and VSP (Fig. 5). Although all of the emergent patterns performed realistically in a qualitative way (i.e. no unusual model behavior or obviously unrealistic patterns across the entire parameter space tested), due to poor confidence in observed data for a number of patterns, we limited our quantification of model robustness to the following patterns: the probability of alpha females reproducing, probability of natal betas, yearling pack size, and beta pack size (in order of priority and level of confidence). We generated additional mean values at a finer scale (VSP and MSP : 0 ? 1 by 0.01) for these four patterns using the regression models described above. We then calculated the percentage of regression estimated means that fell within one standard deviation of mean observed values (Table 2). Our most reliable estimate, the probability of an alpha female denning, was also the most robust to nuisance parameter uncertainty with complete support for all VSP and MSP parameter combinations (i.e. for Random: 100% and Priority: 100% for Prepro = 0.75-0.99). Robustness to deviations in nuisance parameters decreased with decreasing confidence in our observed patterns: the probability of natal betas (Random: 52% and Priority: 48% for Pbeta = 0.15-0.49), yearling pack size (Random: 48% and Priority: 43% for mean Nyearling = 0.50?0.93), and beta pack size (Random: 30% and Priority: 23% for mean Nbeta = 0.46-0.75). Yet, we know pack size is likely biased low due to the use of minimum counts from the red wolf data, especially with respect to the number of betas. Thus, since the model predicted slightly larger numbers of betas and yearlings than 'observed', model robustness was likely stronger than is captured here (i.e. 10% and 8% increase in robustness for pack composition patterns and simultaneous fit respectively with a 0.25 increase in beta and yearling pack size). Since the emergent properties in the IBM did deviate from truth 42 with VSP and MSP, we calibrated the model by selecting values of VSP and MSP that produced model results that matched observed data from the red wolf reintroduction before proceeding with model simulations (see 3.2). 3.2 Efficacy of management in maintaining a stable red wolf population to 200 years Since we were specifically interested in the effects of non-wolf capture effort on the red wolf population, we evaluated the relative extinction dynamics across a range of management effort (i.e. management capture rates per month and seasonal effort) with zone management as currently implemented (i.e. Zone 1: euthanasia of all non-wolf canids, Zones 2 & 3: sterilization of all non-wolf canids). In an effort to reduce the parameter space being tested, as well as to best reflect the system dynamics of the red wolf population, we calibrated nuisance parameters using observed red wolf pack dynamics and pack-level emergent properties (in order to avoid potential problems with spatially-dependent properties) under random (VSP = 0.50, MSP = 0.13, and VMP = 0.50) and priority mate replacement procedures (VSP = 0.45, MSP = 0.13, and VMP = 0.50). Also known as pattern-oriented inverse modeling (Grimm and Railsback 2005, Kramer- Schadt et al. 2007), the above calibration effort selected the most likely parameter combination by filtering out 'unrealistic' parameter combinations by simultaneously fitting multiple emergent patterns to observed system processes (see last paragraph in 3.1 ; performed at the same time as the robust analysis). Since we know management as implemented in the model would have minimal impact on model predictions with full assortative mate replacement (and we know full assortative is not occurring with the actual population), a full assortative procedure was not evaluated in a management context. Only emergent, population-level model output important in assessing persistence potential were collected (i.e. final population size, time to extinction, and probability of extinction) and summarized from the spatial core of the model (i.e. 30 'eastern' 43 packs). Simulations consisted of ?200 iterations for each combination of seasonal effort and capture probability, with nuisance parameters held constant at calibrated values. Since the model predicted a mean time to extinction beyond 50 years in the absence of management, we chose a longer projection length of 200 years in an effort to better understand the influence of management on population persistence. Model predictions were remarkably similar between random and priority mate replacement procedures due to calibration efforts and likely the size of the core assessment area. Thus, only the relative population dynamics of models with red wolf priority mate replacement procedure will be discussed here. In general, predictions of mean final red wolf population size at 200 years increased non-linearly with increases in the probability of monthly non-wolf capture (Fig. 7). Seasonal effort (single season to year round) influenced the rate at which red wolf population size stabilized with increases in monthly capture probability (i.e. populations with more seasonal effort asymptote at lower capture probabilities). However, the sizes at which populations would asymptote were generally consistent across seasonal effort with an average between 70 - 80 individuals (~15 - 16 packs) and was reached by capture probabilities of 0.50 in all cases (Fig. 7). A comparison of single season management effort showed slightly lower predictions of final population size for spring trapping (April ? June) relative to the other three seasons at capture probabilities > 0.25 (Fig. 7a). All other seasons within single season management models produced similar predictions of a final population size. The differences between multi-season models was largely visible between capture probabilities of 0.00 ? 0.30, where increased seasonal effort sustained significantly larger red wolf populations (Fig. 7b). Relative extinction probabilities across management capture probabilities were strongly correlated with mean final population size of red wolves with larger mean populations 44 corresponding to more robust or viable populations. In general, all populations went extinct by year 200 (~39% at 50 years) in the absence of management (Fig. 7; MCP = 0.00). Additionally, nearly complete extinction of all populations occurred at or below capture probabilities of 0.10, 0.05, and 0.00 for single season, 2-season, and 3- and 4-season models respectively. As the probability of capture increased under management, the likelihood of extinction at 200 years decreased until near 0% at a monthly capture probability of approximately 0.30 with minimal seasonal effort (i.e. single season). Seasonal effort also increased the mean time to extinction and influenced the rate at which extinction probabilities approached zero with incremental increases in coyote/hybrid capture probability (i.e. the more trap seasons the lower probability of monthly capture required to attain <0.05% of extinction at any time scale). 3.3 Efficacy of Sterilization vs. Euthanasia Following a similar protocol as in section 3.2, we evaluated red wolf persistence potential (i.e. final population size, probability of extinction, and time to extinction) in the spatial core of the model after varying the 3-zone management theme from its current implementation to euthanasia of non-wolf individuals in all zones. The goal of this modeling effort was to determine the most appropriate zone-specific management based on relative emergent model dynamics. Henceforth, Zone 1, Zones 1-2, and Zones 1-2-3 correspond to the zones for which all non-wolves captured were euthanized. All captures in the remaining zones resulted in sterilization. To improve computation time, we used a coarser parameter space for management effort (monthly capture and seasonal effort) and ran ?200 iterations of 200 year projections for each management combination. Although zone management was likely most subject to variation in spatial scale, space was kept constant as described in section 2.4.6. Simulation results suggested zone management significantly influenced predictions of 45 mean final population size at 200 years (Fig. 8). At moderate to high monthly capture probabilities, populations stabilized around approximately 75 (15 packs), 95 (19 packs), and 150 (30 packs) red wolf individuals under Zone 1, Zones 1-2, and Zones 1-2-3 management respectively. However, at very low monthly capture probabilities as populations climbed out of complete extinction, zone management with increased sterilization effort predicted slightly larger populations on average (Fig. 8). Thus, based on mean final population size, preference for Zone 1 management rapidly switched to Zones 1-2, and in turn Zones 1-2-3, across a seasonal effort dependent range of capture probabilities (Single-season effort : 0-0.20, Two-season effort : 0- 0.15, and 4-season effort : 0-0.10), with increased euthanasia of non-wolf individuals favored at moderate to high monthly capture probabilities (>0.20). We evaluated extinction risk at 50 as well as 200 years because of greater precision in the former (less error propagation). However, we were also interested in how management influenced the long-term stability of the red wolf population, requiring the longer projects of the latter. Extinction probabilities at 50 years under the three zone management scenarios tested were not significantly different from one another after accounting for capture probability and seasonal effort (results using logistic regression with zone, capture probability and season as predictors and binary extinction as a response), suggesting zone had little effect on the risk of extinction at 50 years. However, the model predicted slightly lower probabilities of extinction at 200 years for Zone 1 management (as currently implemented within the recovery area) relative to the other zone management scenarios at monthly capture probabilities less than a seasonal effort dependent threshold (Fig. 9). Once capture probabilities passed beyond this threshold, increased euthanasia (i.e. Zones 1-2 and Zones 1-2-3 management) became equally or more effective at reducing extinction risk until capture probabilities were sufficiently high to stabilize extinction 46 risk at zero for all three scenarios. Finally, for those populations that went extinct, there was no effect of zone management on mean time to extinction at 50 or 200 years. 3.4 Impacts of management on hybridization rates Following a similar protocol as in section 3.3, we evaluated the efficacy of management effort (i.e. capture, zone, and seasonal effort) at reducing the rate of new pairings between wolves and non-wolf canids. However, due to the quantity of data produced in this simulation, we limited our projection length to 50 rather than 200 years. In addition, we tracked the number of red wolf and hybrid individuals (as well as packs) outside the 'recovery' area to assess the influence of management effort on population expansion in red wolves. Model performance suggested monthly non-wolf capture probability had the strongest effect on the annual rate of new hybrid pairings (at 50 years), henceforth ?hybridization rate?. The annual hybridization rate averaged around 11 new hybrid pairings in the absence of management (i.e. capture probability = 0.0). In general, monthly capture probability had a negative non-linear effect with substantial declines in hybridization rates at low capture probabilities and more moderate declines to near stable rates at higher capture probabilities(Fig. 10). Zone management largely influenced the 'stable' point for hybridization rates at moderate to high capture probabilities: Zone 1 (~2.06-2.32), Zones 1-2 (~1.58-1.61), and Zones 1-2-3 (~0- 1.29). Seasonal effort had a weaker, but important effect on hybridization rates. In general, seasonal effort influenced the rate of decline in hybridization with increases in capture probability; as shown in figure 10, more effort equated to a more rapid decline in hybridization rates with increases in capture capture probabilities. However, seasonal effort had an added interaction with Zones 1-2-3 management (i.e. all euthanasia) whereby reduced effort slightly increased the rate of hybridization at low capture probabilities (Fig. 10). This spike in 47 hybridization rate gradually disappeared with increased seasonal effort. Variation in management effort appeared to have minimal impact on red wolf dispersal outside the recovery area, with the model predicting fewer than five red wolves on average within the buffer zone at 50 years (Table 3). However, management did influence the number and proportion of hybrid individuals within the buffer zone (Table 3). The model predicted non-linear declines in hybrid proportions with increases in monthly non-wolf capture probability, eventually stabilizing at moderate to high capture probabilities. The zone scenarios influenced the initial rate of decline in hybrid proportions, as well as the stable point, with increased euthanasia effort favoring slower declines and higher stable hybrid proportions. Increased seasonal effort served to reduce hybrid proportions in the buffer zone between MCPs of 0.00 ? 0.10, but only weakly so. 3.5 Elasticity analysis The elasticities of viability statistics (population size and extinction probability) and hybridization rates to changes in model parameters were assessed using manual perturbation (Morris and Doak 2002) in an effort to identify parameters and/or vital rates potentially important in maintaining red wolf population persistence and reducing the deleterious effects of hybridization. The perturbation analysis was conducted with only Zone 1 management (current recovery area implementation) at no (0.0), moderate (0.5), and high (1.0) monthly capture probabilities crossed with fall, spring and fall, and year-round effort. Individual vital rates were perturbed by 5% while all remaining vital rates were kept constant at predefined values. Since extinctions only occurred at low monthly capture probabilities, the elasticities for the probability of extinction reflect the relative importance of vital rates with minimal capture effort (i.e. elasticities were 0 for all vital rates at moderate and high monthly capture regardless of seasonal effort). The elasticities suggested that alpha wolf mortality, followed by the mean 48 number of wolf recruits per female, had the biggest impact on extinction risk. However, mean final population size and mean final number of wolf packs was most affected by changes in the mean number of wolf recruits, followed by alpha wolf mortality. Increased seasonal effort had the interesting effect of increasing the importance of adult coyote and beta wolf mortality relative to alpha wolf mortality and wolf reproductive rates. Finally, the elasticities for hybridization rate suggested that the annual number of new hybrid pairings was particularly sensitive to adult wolf survival at low capture probabilities and seasonal effort. Whereas at high capture effort, coyote pup mortality became the most influential vital rate in determining annual hybridization rates. 4 Discussion In 1999, the results from a preliminary red wolf PVA suggested the biggest threat to red wolf recovery was hybridization with coyotes. In an effort to evaluate the efficacy of the coyote/hybrid management that followed, our individual-based model was used to evaluate the influence of management on population persistence using models with random and red wolf priority mate replacement procedures calibrated to observed population dynamics. Although mate replacement procedures were coded very differently (i.e. draw from different pools of individuals), model performance was comparable due to efforts to calibrate the model using observed data. In general, the interpretation of model performance within a management framework remained consistent regardless of mate replacement procedure. Assuming the dynamics of the red wolf population were accurately modeled herein, the IBM could provide managers with benchmarks for the minimum proportion of coyotes captured within the recovery area per unit time needed to prevent a decline in red wolf numbers. Model predictions suggest a minimum annual target between 71.8 ? 72.5% (Monthly captures: All year = 0.10 and Single Season = 0.35) in order to maintain a stable red wolf population. 49 We also evaluated alternative zone management strategy by progressively increasing euthanasia effort through zone 3. Model performance suggested zone management as currently implemented with sterilization rather than euthanasia in zones 2 and 3 is more effective at reducing the probability of extinction with single-seasonal effort and low monthly non-wolf capture probabilities. Essentially, when capture likelihood was low and trapping efforts occurred during only a single season each year, the model favored a scenario where an animal was captured, sterilized, and replaced over complete removal thus, preventing a reproductively intact non-wolf (with a low probability of capture) from pair bonding. Interestingly, this effect was in essence the intent of the USFWS by maintaining sterile coyotes/hybrids as non-reproductive placeholders for later displacement by red wolves. However, in our model if capture probabilities were sufficiently high, using euthanasia in zones 2 and 3 was as or more effective than sterilization in reducing probability of extinction, while maintaining larger red wolf population numbers and reducing hybridization rates (Figs. 8 and 10). Intuition suggests that in a red wolf population managed for expansion, surplus individuals would disperse outside the recovery area and interbreed with coyotes. However, the model predicts very little establishment of red wolves beyond the recovery area with zone management as currently implemented (regardless of effort) and a landscape saturated with coyotes (<5 red wolf individuals on average at 50 years). Supporting our model predictions, Bohling et al. (in press) found little evidence of red wolf establishment outside the recovery area. However, Bohling et al. also found little evidence of red wolf genetic introgression in the coyote population immediately west of areas managed for red wolves. Although genetics were not modeled explicitly in the IBM, the model contradicts the latter findings of Bohling et al. and supports the idea that some introgression of red wolf genetics into the coyote population 50 immediately west of the recovery area would occur under management, at least on average (Table 3). These findings pose the question then as to why we do not see more movement of red wolf genetics outside the recovery area. Could management within the recovery area be reducing the rate of hybridization outside the recovery area or is it simply chance that has limited introgression of red wolf genetics to date? Either way, these findings highlight the need to further understand red wolf dispersal and density-dependent demographics, specifically mate replacement processes between and within these two canid species. In an effort to better understand red wolf viability with coyotes on the landscape, a viability model for this two-species system necessitated greater realism and more informed red wolf demographics. One drawback to building a complex model is the difficulty in parameterizing such a model. This problem was mediated in our case owing to 20 years of extensive data collection on the red wolf reintroduction. Nonetheless, it was necessary to include three nuisance parameters to fill gaps in our understanding of the system, while maintaining the most biologically appropriate and parsimonious model structure possible. Vagrant movement probability (VMP) was the only one of three nuisance parameters that did not influence emergent system dynamics in a biologically significant way. The degree to which the remaining two, vagrant settlement probability (VSP) and mate selection probability (MSP), influenced model predictions demonstrates the need to better understand mate replacement and settlement processes within coyotes and red wolves. Although we are confident our findings accurately reflect the relative dynamics of the system following model calibration, they do not preclude the need for more informed coyote demographics from the region, as well as specifically a better understanding of density-dependent mate replacement processes in both canids. Overall, perturbations in the model's nuisance parameters produced expected, as well as 51 biologically appropriate, outcomes. However, the apparent interaction between vagrant settlement and mate replacement and the interaction's influence on the probability of an alpha female reproducing was one surprising emergent property of the model. In particular, there were substantial declines in the probability of denning at high vagrant settlement and low mate selection probabilities (Fig. 6). Although somewhat counterintuitive, high vagrant settlement within established packs effectively reduced the pool of available mates. The interaction suggests that highly selective individuals (low MSP) may substantially reduce reproductive rates within a population when vagrant settlement probabilities are high. This result was similar across all three mate replacement procedures, but was most notable with random mate replacement where the mate pool experienced the greatest change in size from low to high settlement probability. Although interesting, it is reasonable to assume that within most real systems, mate selectivity changes with density (i.e. individuals become less selective at lower densities); a phenomenon that we did not model. Another interesting result was that increases in the probability of a vagrant settling within an established pack (VSP) had an indirect negative effect on new pack establishment (Fig. 3). Opportunities for vagrant individuals to meet and establish new packs was reduced due to a smaller vagrant pool associated with higher vagrant settlement probabilities, and in part contributed to the heteroschedasticity in final population size shown in figure 4. These predictions suggest that populations undergoing high rates of settlement in established social groups could impede population expansion, especially within a competitive two species framework such as the red wolf and coyote system where selection may favor the species with a higher yield of surplus individuals who in turn can quickly occupy available space. In the same context, actual red wolf data suggests a high degree of unrelatedness amongst known pack 52 members (Mahoney unpublished results). In our model, the probability of betas belonging to their natal pack, loosely analogous to intrapack relatedness, was one of the emergent properties used to calibrate the model with a moderate vagrant settlement rate of VSP ? 0.45. In turn, wolves within socially disruptive systems (i.e. high rates of alpha and/or beta mortality) have been shown to experience a higher degree of unrelatedness amongst pack members relative to undisturbed systems (Rutledge et al. 2010b). With respect to management, perhaps encouraging pack retention of related individuals among non-alpha pack members, whether or not by reducing socially disruptive processes in the recovery area, may discourage canid settlement within existing packs and encourage new pack establishment in available or coyote-occupied territories. In conclusion, the current IBM provides some important information regarding the management of the reintroduced red wolf population, as well as some interesting insights into factors affecting population viability of species within competitive systems in general. The IBM provides the foundation for future red wolf viability modeling and could, for example, be further enhanced by incorporating space use behavior, allowing us to evaluate secondary red wolf release sites. In addition, the model could also be used to explore competitive interference disequilibrium between red wolves and coyotes in an open system (as opposed to the 'closed' system modeled here). Finally, individual genotypes could be included in the model, facilitating the exploration of questions pertaining to heredity and introgression. Ultimately, the model will hopefully evolve to assist the efforts of management officials towards more informed management and to help ensure the continued success of red wolves in the wild. 53 Table 1: An outline for the IBM detailing sub-model structure, vital rates, and stochastic processes. aStochastic process describes the method by which stochasticity was incorporated in the model. Demographic stochasticity was the most common form and used uniform distributions with derived means in a Bernouilli Trial. Initial beta pack size and number of recruits accounted for some environmental variation during model fitting by using all 20 years of red wolf data. bNvag ? Number of vagrants in population (both), P ? Number of packs within population (both), A ? Age in months. clogit equation used in a logistic model for determining dispersal probability : Npack ? Pack size, TDis ? Time in months since last dispersal, D ? Dominant status (logical positive), NDis ? Number of prior dispersals, SFallWin ? Fall or winter (logical positive). 54 Table 1 : Description Value/Equation Model Initialization Number of territories in core assessment area 30 (5x6) Number of territories in buffer zone 30 (5x6) Beta pack size Sex ratio 0.5000 Vital Rates ? Red Wolf Mortality in alphas 0.0072 Uniform[0,1] Mortality in betas 0.0295 Uniform[0,1] Mortality in yearlings 0.0170 Uniform[0,1] Mortality in pups 0.0168 Uniform[0,1] Mortality in vagrants 0.1178 Uniform[0,1] Recruits per alpha female Mortality in adults Uniform[0,1] Mortality in pups 0.1 Uniform[0,1] Uniform[0,1] Pups per alpha female Social Transition Alpha replacement (Mate Selection, MSP) 0.1 Uniform[0,1] Uniform[0,1] Vagrant movement (VMP) 0.0 ? 1.0 Vagrant settlement (VSP) 0.0 ? 1.0 Management (core assessment area only) Monthly non-wolf capture probability (MCP) 0.0 ? 1.0 Capture seasons single season through year-round Stochastic Processa ?=1.197 P0=0.425 Zero-inflated Poisson (fitted to data from 1987 ? 2007) ?=2.375 k=8.579 P0=0.248 Zero-inflated Negative Binomial (fitted to data from 1987 ? 2007) Vital Rates ? Coyotes (Pitt et al. 2003,Conner et al. 2008) 0.01000-0.0003A+0.00025A2 Mortality in vagrantsb [0.008+0.089(Nvag/P)]-0.0003A+0.00025A2 6.93-0.72Npack Dispersal ? Red Wolvesc logit(-5.0967-0.0382A-0.0120TDis-1.2193D+1.5850N Dis+0.9317SFallWin+0.8258Npack-0.0692Npack 2) Dispersal ? Coyotes (Pitt et al. 2003,Conner et al. 2008) 0.005Npack 2 Table 2: Results from the robust analysis showing the percentage of parameter combinations capturing observed system patterns (within 1 SD of mean). The robustness of yearling and beta pack size, and in turn simultaneous fit, are likely biased low due to the use of minimum counts for each. 55 Table 2 : Single pattern fit (%) Random 10201 100 52 48 30 2 (n=236) Red Wolf Priority 10201 100 48 43 23 1 (n=145) Parameter Combinations Simultaneous fit (%) Mate Selection Procedure Probability of Denning Natal beta proportion Yearling Pack Size Beta Pack Size Table 3: Estimated canid proportions inside and outside the core assessment area with Spring and Fall management captures and under priority mate selection. MCP : Management non-wolf capture probability. 56 Table 3 : Proportion of Hybrids Proportion of Red Wolves inside core area inside core area MCP=0 MCP=0.5 MCP=1 MCP=0 MCP=0.5 MCP=1 Zone 1 0.81 0.14 0.01 0.16 0.72 0.97 Zones 1-2 0.78 0.12 0.01 0.21 0.79 0.98 Zones 1-2-3 0.71 0.09 0 0.33 0.87 0.99 outside core area outside core area MCP=0 MCP=0.5 MCP=1 MCP=0 MCP=0.5 MCP=1 Zone 1 0.82 0.72 0.60 0.001 0.002 0.005 Zones 1-2 0.84 0.75 0.64 0.002 0.003 0.007 Zones 1-2-3 0.98 0.97 0.94 0.006 0.013 0.029 Euthanasia Zones Figure 1: A diagrammatic representation of model structure and scheduling within the red wolf and coyote individual-based model. ARNWR : Alligator River National Wildlife Refuge. 57 Figure 2: Surfaces depicting predicted means from negative binomial (A-C) and logistic (D) regression models fitted to emergent pack compositional characteristics: (A) pack size, (B) number of yearlings within a pack, (C) number of betas within a pack, and (D) probability of a beta belonging to it's natal pack (D). In all cases, the red wolf priority mate replacement procedure was used. 58 Figure 3: A surface depicting the predicted mean number of vagrants from a negative binomial regression model fitted to vagrant settle probabilities and mate selection probabilities with the red wolf priority mate replacement procedure. 59 Figure 4: Distributions of model output for red wolf population size at 50 years for various levels of VSP (Vagrant Settlement Probability) and MSP (Mate Selection Probability). In all cases, the red wolf priority mate replacement procedure was used. 60 Figure 5: Surfaces depicting the predicted annual mean number of new hybrid packs from a negative binomial regression model under (A) random and (B) red wolf priority mate replacement procedures. 61 Figure 6: Surface depicting the predicted probability of an alpha female reproducing from a logistic regression model with the red wolf priority mate replacement procedure. 62 Figure 7: A comparison of the effects of capture probability and seasonal effort on the mean final red wolf population size at 200 years with the red wolf priority mate replacement procedure (Calibrated Values: VSP = 0.45 and MSP = 0.13). Error bars represent 65% rank intervals. 63 Figure 8: Mean final red wolf population size at three levels of zone management and four seasonal efforts with the red wolf priority mate replacement procedure. Error bars represent 65% rank intervals around means at 200 years. 64 Figure 9: Probability of extinction (200yrs) at three levels of zone management and four seasonal efforts with the red wolf priority mate replacement procedure. 65 Figure 10: Mean number of new hybrid pairings annually at three levels of zone management and four seasonal efforts with the red wolf priority mate replacement procedure (Calibrated Values: VSP = 0.45 and MSP = 0.13). Error bars represent 65% rank intervals around means at 50 years. 66 References Adams, J.R., C. Lucash, L. Schutte, and L.P. Waits. 2007. Locating hybrid individuals in the red wolf (Canis rufus) experimental population area using a spatially targeted sampling strategy and faecal DNA genotyping. Molecular ecology 16:1823-1834. Ak?akaya, H.R. 2002. RAMAS GIS: linking spatial data with population viability analysis (version4.0). Applied Biomathematics. Setauket, N.Y. Beissinger, S., and D.R. McCullough, eds. 2002. Population Viability Analysis. University of Chicago Press, Chicago, Illinois, USA. Beissinger, S., and M. Westphal. 1998. On the use of demographic models of population viability in endangered species management. Journal of Wildlife Management 62:821-841. Bohling, J.H., and L.P. Waits. unpublished. Assessing the prevalence of hybridization between sympatric Canis species in a region surrounding the red wolf (Canis rufus) recovery area in North Carolina. Boyce, M. S. 1992. Population Viability Anlaysis. Annual Review of Ecology and Systematics 23:481-506. Brook, B.W. 2000. Pessimistic and optimistic bias in population viability analysis. Conservation Biology 14:564-566. Brook, B.W., J.J. O'Grady, A.P. Chapman, M.A. Burgman, H.R. Akcakaya, and R. Frankham. 2000. Predictive accuracy of population viability analysis in conservation biology. Nature 404:385-387. Brook, B.W., M.A. Burgman, and H.R. Akcakaya. 2002. Critiques of PVA ask the wrong questions: Throwing the heuristic baby out with the numerical bath water. Conservation Biology 16:262-263. Burgman, M. and H. Possingham. 2000. Population Viability analysis for conservation: the good, the bad and the undescribed. In: Young, A.G. and G.M. Clarke (eds). Genetics, demography and viability of fragmented populations. Cambridge University Press pgs. 97-112. Caswell, H. 2001. Matrix Population Models: Construction, Analysis, and Interpretation, second ed. Sinauer Associates, Massachusetts. Caughley, G. 1994. Directions in conservation biology. Journal of Animal Ecology 63:215-244. 67 Conner, M. M., M. R. Ebinger, and F. F. Knowlton. 2008. Evaluating coyote management strategies using a spatially explicit, individual-based, socially structured population model. Ecological Modelling 219:234-247. Coulson, T., G.M. Mace, E. Hudson, and H. Possingham. 2001. The use and abuse of population viability analysis. Trends in Ecology & Evolution 16:219-221. Doak, D.F., K. Gross, and W.F. Morris. 2005. Understanding and predicting the effects of sparse data on demographic analyses. Ecology 86:1154-1163. Fredrickson, R. J., and P. W. Hedrick. 2006. Dynamics of hybridization and introgression in red wolves and coyotes. Conservation Biology 20:1272-1283. Gaona, P., P. Ferreras, and M. Delibes. 1998. Dynamics and viability of a metapopulation of the endangered Iberian lynx (Lynx pardinus). Ecological Monographs 68:349-370. Grimm, V., U. Berger, F. Bastiansen, S. Eliassen, V. Ginot, J. Giske, J. Goss-Custard, T. Grand, S. K. Heinz, G. Huse, A. Huth, J. U. Jepsen, C. Jorgensen, W. M. Mooij, B. Muller, G. Pe'er, C. Piou, S. F. Railsback, A. M. Robbins, M. M. Robbins, E. Rossmanith, N. Ruger, E. Strand, S. Souissi, R. A. Stillman, R. Vabo, U. Visser, and D. L. DeAngelis. 2006. A standard protocol for describing individual-based and agent-based models. Ecological Modelling 198:115-126. Grimm, V., U. Berger, D.L. DeAngelis, J.G. Polhill, J. Giske, and S.F. Railsback. 2010. The ODD protocol: A review and first update. Ecological Modelling 221:2760-2768. Grimm, V., N. Dorndorf, F. Frey-Roos, C. Wissel, T. Wyszomirski, and W. Arnold. 2003. Modelling the role of social behavior in the persistence of the alpine marmot Marinota marmota. Oikos 102:124-136. Grimm, V. and S.F. Railsback. 2005. Individual-based Modeling and Ecology. Princeton University Press, Princeton, New Jersey. Grimm, V., and C. Wissel. 2004. The intrinsic mean time to extinction: a unifying approach to analysing persistence and viability of populations. Oikos 105:501-511. Gusset, M., O. Jakoby, M. S. Muller, M. J. Somers, R. Slotow, and V. Grimm. 2009. Dogs on the catwalk: Modelling re-introduction and translocation of endangered wild dogs in South Africa. Biological Conservation 142:2774-2781. Kays R., A. Curtis, and J.J. Kirchman. 2010. Rapid adaptive evolution of northeastern coyotes via hybridization with wolves. Biol Letters 6:89?93. 68 Kelly, B.T., P.S. Miller, and U.S. Seal, eds. 1999. Population and habitat viability assessment workshop for the red wolf (Canis rufus). Conservation Breeding Specialist Group (CBSG, SSC/IUCN). Kelly, B. T. 2000. Red wolf recovery program adaptive work plan - FY00?FY02. United States Fish and Wildlife Service Report, Manteo, North Carolina, USA. Knowlton, F.F., E.M. Gese, J.R. Adams, K. Beck, T.K. Fuller, D.L. Murray, T. Steury, M.K. Stoskopf, W. Waddell, and L.P. Waits. unpublished. Managing hybridization in endangered species recovery: The case of the red wolf. Kramer-Schadt, S., E. Revilla, T. Wiegand, and U. Breitenmoser. 2004. Fragmented landscapes, road mortality and patch connectivity: modelling influences on the dispersal of Eurasian lynx. Journal of Applied Ecology 41:711-723. Kramer-Schadt, S., E. Revilla, T. Wiegand, and V. Grimm. 2007. Patterns for parameters in simulation models. Ecological Modelling 204:553-556. Lacy, R.C. 2000. Considering threats to the viability of small populations. Ecological Bulletins 48:39-51. Mech, L.D., and L. Boitani (Editors). 2003. Wolves: Behavior, Ecology, and Conservation. University of Chicago Press, Chicago, IL. Miller, P.S., and R.C. Lacy. 2005. VORTEX. A stochastic simulation of the simulation process. Version 9.50 user's manual. Conservation Breeding Specialist Group (IUCN/SSC). Apple Valley, Minnesota. Morris, W. F., and D.F. Doak. 2002. Quantitative Conservation Biology: Theory and Practice of Population Viability Analysis. Sinauer Associates, Massechussetts. Nowak, R. M. 2002. The original status of wolves in eastern North America. Southeastern Naturalist 1: 95-130. Phillips, M.K., V.G. Henry, and B.T. Kelly. 2003. Restoration of the red wolf. Pp. 272-288 In Mech, L.D., and L. Boitani, eds. Wolves: Behavior, ecology, and conservation. University of Chicago Press, Chicago, IL. Pitt, W.C., P.W. Box, and F.F. Knowlton. 2003. An individual-based model of canid populations: modelling territoriality and social structure. Ecological Modelling 166:109-121. Possingham, H.P., D.B. Lindenmayer, and G. Tuck. 2002. Decision theory thinking for 69 population viability analysis. In: Beissinger, S.R. And D.R. McCullough (eds). Population viability analysis. University of Chicago Press. Rutledge, L.Y., B.R. Patterson, K.J. Mills, K.M. Loveless, D.L. Murray, and B.N. White. 2010a. Protection from harvesting restores the natural social structure of eastern wolf packs. Biological Conservation 143:332-339. Rutledge, L.Y., C.J. Garroway, K.M. Loveless, and B.R. Patterson. 2010b. Genetic differentiation of eastern wolves in Algonquin Park despite bridging gene flow between coyotes and grey wolves. Heredity 1-12. Roth, J.D., D.L. Murray, and T.D. Steury. 2008. Spatial dynamics of sympatric canids: modeling the impact of coyotes on red wolf recovery. Ecological Modelling 214:391-403. Swanack, T. M., W.E. Grant, and M.R. J. Forstner. 2009. Projecting population trends of endangered amphibian species in the face of uncertainty: A pattern-oriented approach. Ecological Modelling 220:148-159. Steury, T.D., K. Beck, A. Byers, and D.L. Murray. In review. Does density-dependent habitat use limit the utility of resource selection: models? Ecological Applications. Thiele JC and V. Grimm. 2010. NetLogo meets R: Linking agent-based models with a toolbox for their analysis. Environmental Modelling and Software 25: 972-974. U.S. Fish and Wildlife Service. 1989. Red Wolf Recovery Plan. U.S. Fish and Wildlife Service, Atlanta, Georgia. U.S. Fish and Wildlife Service. 1990. Red Wolf Recovery/Species Survival Plan. U.S. Fish and Wildlife Service, Atlanta, Georgia. U.S. Fish and Wildlife Service. 2007. Red wolf 5-year status review: summary and evaluation. U.S. Fish and Wildlife Service, Atlanta, Georgia. Vucetich, J. A., and S. Creel. 1999. Ecological interactions, social organization, and extinction risk in African wild dogs. Conservation Biology 13:1172-1182. Wakamiya, S. M., and C.L. Roy. 2009. Use of monitoring data and population viability analysis to inform reintroduction decisions: Peregrine falcons in the Midwestern United States. Biological Conservation 142:1767-1776. Wiegand, T., J. Naves, T. Stephan, and A. Fernandez. 1998. Assessing the risk of extinction for the brown bear (Ursus arctos) in the Cordillera Cantabrica, Spain. Ecological 70 Monographs 68:539-571. Wiegand, T., F. Jeltsch, I. Hanski, and V. Grimm. 2003. Using pattern-oriented modeling for revealing hidden information: a key for reconciling ecological theory and application. Oikos 100:209-222. Wilensky, U. 1999. NetLogo. http://ccl.northwestern.edu/netlogo/. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL. 71