Variation in vital-rate sensitivity between populations of Texas horned lizards

Demographic studies of imperiled populations can aid managers in planning conservation actions. However, applicability of findings for a single population across a species’ range is sometimes questionable. We conducted long-term studies (8 and 9 years, respectively) of 2 populations of the lizard Phrynosoma cornutum separated by 1000 km within the historical distribution of the species. The sites were a 15-ha urban wildlife reserve on Tinker Air Force Base (TAFB) in central Oklahoma and a 6000-ha wildland site in southern Texas, the Chaparral Wildlife Management Area (CWMA). We predicted a trade-off between the effect of adult survival and fecundity on population growth rate (λ), leading to population-specific contributions of individual vital rates to λ and individualized strategies for conservation and management of this taxon. The CWMA population had lower adult survival and higher fecundity than TAFB. As predicted, there was a trade-off in the effects of adult survival and fecundity on λ between the two sites; fecundity affected λ more at CWMA than at TAFB. However, adult survival had the smallest effect on λ in both populations. We found that recruitment in P. cornutum most affected λ at both sites, with hatchling survival having the strongest influence on λ. Management strategies focusing on hatchling survival would strongly benefit both populations. As a consequence, within the constraint of the need to more accurately estimate hatchling survival, managers across the range of species such as P. cornutum could adopt similar management priorities with respect to stage classes, despite intra-population differences in population vital rates.


Introduction
Effective strategies for management and recovery of imperiled species are a central focus of conservation biology. Knowledge of interspecific variation in vital rates (and associated life histories) can help develop generalized management and conservation priorities for species lacking solid demographic information (Saether et al. 1996;Heppell et al. 2000). Demographic modeling is an important part of this process, with perturbation analyses serving as a standard tool (Caswell 2000). Prospective perturbation analyses (Caswell 2000), such as sensitivity and elasticity, examine the potential effects of changes in a vital rate (e.g., adult survival) on a dependent variable [e.g., population growth rate (k)]. Life-stage simulation analysis (LSA), a form of retrospective perturbation analysis (Wisdom et al. 2000), complements sensitivity/elasticity analyses by including variation (either stochasticity in the system or uncertainty in vital rate estimates) in perturbation analyses. A combination of these analyses is valuable for planning species conservation and recovery efforts, as they help managers decide which aspects of a species' life history to target for conservation action, improving the effectiveness of these actions.
Intraspecific variation in demographic parameters has received little attention (Fredericksen et al. 2005) relative to interspecific and within-population studies of life histories. However, such variation leads to associated differences in the elasticity patterns of those rates (Oli and Dobson 2003), with concomitant differences in their relative importance to the population growth rate. The result is that conspecific or congeneric populations may have different demographic responses to environmental perturbations and stressors, such as climate change (e.g., Lagopus sp., Sandercock et al. 2005; Martin 2012)o r harvest management (Servanty et al. 2011;Devenish-Nelson et al. 2013). Intraspecific variation in vital rates also has led to recent conclusions that effective management or conservation may require population-or situation-specific strategies (Johnson et al. 2010;Servanty et al. 2011;Devenish-Nelson et al. 2013). For example, Johnson et al. (2010) found that the relative importance of individual vital rates in affecting population growth differed among populations of endangered Ovis canadensis sierra, depending on the magnitude of those vital rates and the population trajectory.
Vital-rate variation among species corresponds to differing life-history strategies, and has long been theorized to represent an evolutionary trade-off between survival and reproduction (Stearns 1976(Stearns , 1977. Lizards have been widely used as models to explore trade-offs among lifehistory traits (Dunham and Miles 1985;Miles and Dunham 1992;Shine and Charnov 1992;Bauwens and Diaz-Uriarte 1997;Clobert et al. 1998). This taxon attracts study because lizards exhibit wide variation in life-history traits such as age at maturity, survival of various life stages, annual and lifetime clutch frequency, and clutch size in relation to body size (Tinkle et al. 1970). More recently, studies of life-history traits of lizards have focused on intraspecific variation to elucidate the connection between such variation and local environmental factors acting on different populations (Sears 2005;Radder 2006;Sperry et al. 2010).
Among squamates, horned lizards (Phrynosoma spp.) are a useful model for examining questions regarding the conservation implications of population variation in vital rates. First, they exhibit considerable variation in life history across their range (Ballinger 1974;Pianka and Parker 1975;Montgomery et al. 2003;Endriss et al. 2007). For example, body sizes among Phrynosoma cornutum vary along a latitudinal gradient; females in central Oklahoma reach maturity at approximately 63 mm snout-vent length (SVL; Endriss et al. 2007); whereas in central Texas, mature females are 70-75 mm (Ballinger 1974). Clutch sizes covary with body size, with a mean clutch size of 17.4 eggs in central Oklahoma (Endriss et al. 2007), compared with 29 eggs in central Texas (Ballinger 1974). Additionally, because many species of horned lizards are endangered, threatened, or in decline (Carpenter et al. 1993;Donaldson et al. 1994;Jennings and Hayes 1994;Grant and Doherty 2007), between-population variation in vital rates have potential consequences for conservation planning.
Assessing which vital rates are most important for restoring vertebrate populations and whether these rates  Sherbrooke (2003) differ between populations is of immediate application in management contexts. P. cornutum are rare in Oklahoma and threatened in Texas, and because they are relatively well-studied, they are an ideal species for exploring these issues with ample opportunity to apply findings to conservation and management. We studied 2 populations of P. cornutum separated by 1000 km along a north-south axis (Fig. 1). Our specific objectives were twofold: (1) to describe vital-rate variation in 2 populations of P. cornutum using long-term demographic data (e.g., 17 total studyyears); and (2) to determine the relative influence of individual vital rates on k in these 2 populations. Based on the known latitudinal variation (Ballinger 1974;Montgomery et al. 2003;Endriss et al. 2007), we predicted the southern population would have lower survivorship and higher fecundity because of the interactions among temperature, activity, and predation risk (Adolph and Porter 1993). In turn, we predicted that this intraspecific variation in vital rates would be manifested in population-specific contributions of individual vital rates to k. Finally, we expected that support of our predictions would lead to populationspecific strategies for conservation and management of this taxon.

Study areas
Our 2 study sites were separated by 1000 km, representing one population in the northern part of the species distribution and another population in the southern part of the distribution (Price 1990). The northern site, Tinker Air Force Base (TAFB; Oklahoma County, OK, USA, 35°24 0 58 00 N, 97°24 0 41 00 W), was a largely urban site on the outskirts of Oklahoma City. Mean annual temperature during 1971-2000 was 15.6°C, with mean January and July temperatures of 2.6 and 27.7°C, respectively. Mean annual precipitation was 91.1 cm for the same years (National Climatic Data Center 2002). Of the 2000-ha Base, approximately 500 ha were natural habitat. These areas were dominated by mixed oak-hardwood forests and a mixture of native and non-native grasslands (see Endriss et al. 2007 for a detailed description of floral community). Research activities were focused on the population of P. cornutum on Wildlife Reserve 3 (WR3) and the surrounding areas. Wildlife Reserve 3 is a natural area (ca. 15 ha) on the southwestern side of TAFB, dominated by grassland with patches of woody vegetation and gravel trails and centered around 2 man-made ponds. As of 2011, WR3 was nearly completely surrounded by residential housing or military buildings. Habitat restoration since 2005 and construction of a military housing development in a 7.4-ha area directly adjacent to WR3 in 2008-2010 have disturbed the Reserve. Management activities designed to restore prairie habitat have included tree removal, disking, mowing, spraying with herbicides, and seeding with native grasses and forbs. Although these activities do not seem to affect P. cornutum space use, there is evidence that they have negatively affected adult survival (Wolf et al. 2013).
The southern site, Chaparral Wildlife Management Area (hereafter CWMA; Dimmit and La Salle counties, Texas, 28°19 0 40 00 N, 99°24 0 39 00 W), was a 6,150-ha area managed by the Texas Parks and Wildlife Department (Hellgren et al. 2010). Mean annual temperature at nearby Catarina, TX during 1971-2000 was 22.1°C, with mean January and July temperatures of 11.9 and 30.5°C, respectively. Mean annual precipitation was 51.8 cm (National Climatic Data Center 2002). The site, which was predominantly honeymesquite woodlands or parklands, experiences managed grazing and prescribed burns (Flanders et al. 2006;Hellgren et al. 2010). A detailed site description is available in Flanders et al. (2006).

Field methods
Field methods were similar between sites unless otherwise noted. Work on WR3 occurred from 2003 to 2011 (Endriss et al. 2007;Wolf et al. 2013), whereas study on CWMA occurred from 1998 to 2005 (Hellgren et al. 2010). We captured lizards during April-August through intensive visual searching, road cruising, and fortuitous encounters, and recorded basic morphometric information for each lizard, including snout-vent length (SVL), total length (TL), mass, and sex. Intensive visual searches consisted of slowly walking back and forth across search areas while looking for lizards. We attempted to evenly and thoroughly search all areas of the WR3 field site except areas where vegetation was so thick and high that we probably would not have been able to detect lizards when they were present. Lizards were hand-captured upon detection. We implanted lizards [5.0 g with a passive integrated transponder (PIT) tag (0.5 g), or clipped a unique combination of toes for smaller lizards. Lizards had a single toe clipped as a secondary mark (Hellgren et al. 2010). We attached a 0.95-1.95 g radio transmitter (BD-2, Holohil Systems Ltd., Carp, Ontario, Canada or SOPR-2038, Wildlife Materials Inc., Murphysboro, IL, USA) to each lizard on WR3 if the transmitter was \10 % of the lizard's body mass. On CWMA, we attached a transmitter-backpack bundle totaling 3 g (\8 % of lizard body mass) to adult lizards (Hellgren et al. 2010). The difference in transmitter mass resulted from larger lizards on CWMA (R. T. Kazmaier, unpublished data). Transmitters were attached by the same methods given in Endriss et al. (2007) and Hellgren et al. (2010), and radio-monitoring followed similar regimes at the prior studies on these sites (Endriss et al. 2007;Hellgren et al. 2010).
We monitored the locations of radiomarked lizards on both sites 1-5 times weekly during the active lizard season (Apr-Aug). On WR3, we also monitored lizards at least biweekly from August until they entered hibernation (generally Oct-Dec). We homed to each lizard's location, which were recorded in Universal Transverse Mercator (UTM) coordinates using the North American Datum 1983 (NAD83). We varied the times during which we tracked lizards each day to obtain a representative sample of locations across all daylight hours. We also observed telemetered lizards to gather data on the subvital rates that compose fecundity (see below). We tracked females 2-3 times daily throughout nesting periods, observed nesting activity, and gathered data on clutch frequency. Because horned lizards spend at least a day digging a nest, we observed and recorded nest locations, and monitored nests daily from the date they were laid to their hatching date. Further details on detecting and monitoring nests are available in Endriss et al. (2007). We defined lizards as hatchlings during the active season in which they hatch, as juveniles between their first and second hibernations, and breeding adults after the second hibernation of a lizard's lifetime.

Model parameter estimates
Life-stage simulation analysis relies on randomly drawing vital-rate values from a probability distribution based on the estimated mean and variance of said vital rate. Ideally, the mean and variance for each vital rate should be estimated from empirical data, and estimates of process variance and sampling variance should be separated. (Morris and Doak 2002). We attempted to use White's (2000) method for variance decomposition, but in all instances, we failed to calculate realistic estimates of process variance, probably because sampling variance was too great relative to process variance (Gould and Nichols 1998). Therefore, we incorporated empirical estimates of variance (process and sampling variance combined) into our LSA wherever possible (Wisdom et al. 2000). We note that including sampling variance will over-estimate the variance of each parameter.
We modeled vital rates limited to a 0-1 scale (survival rates, proportion of females breeding, proportion of females double-clutching, and hatch rates) using b distributions, and modeled clutch sizes using a stretched b distribution (Wisdom et al. 2000;Morris and Doak 2002). The shape of a b distribution is dictated by the parameters a and b; population mean and variance estimates can be transformed to a and b estimates for a b distribution (Morris and Doak 2002). However, if the variance is too great relative to the mean for a parameter (e.g., due to sampling variance, see above), the resulting b distribution often does not match the distribution expected based on biological information for the study species. The distribution may become bimodal, or the distribution may have a mode much closer to 0 or 1 than to the mean. This property of b distributions may be useful in modeling boom-and-bust phenomena, but modelers should be aware of the potential for creating a bimodal distribution when a unimodal distribution is more biologically accurate based on prior knowledge of the study species. This scenario is particularly common when field data yield wide sampling variances, which may or may not reflect process variance accurately. Biologically, these would represent more extreme values or a more skewed distribution than is realistic based on empirical data and known information about the species. This unrealistic skew and/or high frequencies of 0 and 1 occurred in a number of vital rate distributions for P. cornutum at WR3 and CWMA. In several cases, we therefore adjusted the variance for some parameter estimates to better match biological expectations. All means, however, were empirically based.

Fecundity
We defined fecundity as the number of female offspring per female per year. Similar to Wisdom and Mills (1997), we calculated fecundity using a number of subvital rates, with fecundity = sex ratio 9 [(proportion females breeding 9 first clutch size) ? (proportion females doubleclutching 9 second clutch size)] 9 nest survival 9 hatch rate. Following Endriss et al. (2007), we assumed P. cornutum populations to have a proportion of 50 % females. For the purposes of LSA, we fixed the sex ratio instead of drawing a random value from a probability distribution for each model iteration, as this vital rate seemed constant over time (Endriss et al. 2007;A. J. Wolf, unpublished data).
We defined nest survival as the probability of at least one hatchling surviving incubation to hatch. Nests did not survive if the nest was depredated, flooded, or laid incorrectly (i.e., one female dug a nest hole, but laid her eggs outside the hole and back-filled an empty nest, leaving the eggs exposed). Some nests on WR3 had eggs that hatched successfully but whose hatchlings could not escape the nest chamber and subsequently died. We considered these nests to have survived, and incorporated the deaths of hatchlings into the hatch rate (see below). As there was no relationship between hatchling survival (which is often heavily affected by environmental conditions such as temperature and humidity) and laying date (R. T. Kazmaier, unpublished data), we assumed that nest survival rates were the same for first and second clutches. We calculated a weighted mean (White 2000) estimate of nest survival at WR3 using estimates and sample sizes from field data collected there during -2005Endriss et al. 2007Endriss et al. ), 2010Endriss et al. , and 2011. We used the sample variance, which was calculated as the average variance among years. Nest survival at CWMA was estimated based on nests monitored between 1998 and 2001, pooling across years. We calculated sampling variance of this estimate as s 9 (1-s)/n, where s is the estimated probability of nest survival and n is the number of nests in the sample.
We estimated the proportion of adult females who were reproductively active using data from the center of the P. cornutum distribution (Fig. 2 in Ballinger 1974). The proportion we calculated, 0.9636, seemed appropriate for WR3, as all adult females that were telemetered during the first nesting periods of 2010 and 2011 displayed nesting behavior. Because of low survival at CWMA, few females were monitored throughout the long nesting period that occurs in southern Texas. However, the length of the nesting period and the abundance of food (Pogonomyrmex spp.) at the site should promote breeding by nearly all females. We therefore assumed the estimate we obtained from Ballinger (1974) was applicable to CWMA.
The greatest possible variance that could be used in a b distribution with a mean of 0.9636 would be 0.0351. However, random draws from a b distribution with mean = 0.9636 and variance = 0.03 produce a value of 1 almost all the time. Assuming some females may be unable to reproduce in resource-poor years, we used a variance of 0.01 to incorporate a small amount of variance into the model. This combination of mean and variance resulted in values [0.9 for almost 90 % of random draws (n = 10,000).
We determined clutch size after hatching by digging up each nest and counting the number of eggs. We estimated the mean number of eggs in the first clutch laid by females each season using field data from 3 time periods: 2004-2005 combined (Endriss et al. 2007), 2010, and 2011. We used an unweighted mean of the first clutch size for each time period as the estimated mean for the probability distribution and used the variance among the average clutch size of the time periods. We determined the mean clutch size of second clutches (n = 3) during 2011 on WR3 using the methods described for first clutch-size estimates. Because this estimate was based on data from 1 year, we used the variance of the mean estimate from these 3 nests for the LSA. Field methods for estimating clutch size at CWMA were similar to those we used on WR3. However, estimates at CWMA were complicated by continual nest-laying throughout the active season. Therefore, we could not tell which nests were first, second, or even perhaps third clutches. Evidence of triple-clutching in south Texas could not be confirmed at CWMA, so we assumed all clutches were first or second nests. To estimate sizes of first and second clutches, we assumed the total sample mean was a weighted mean of first and second clutches. We used the proportion of females nesting (Ballinger 1974), the proportion double-clutching on CWMA (see below), and the ratio of first: second clutch size on WR3 to solve algebraically for first and second clutch sizes on CWMA. Because variance estimates for clutch size on WR3 were similar between the first and second clutches, we assumed the same for CWMA. We calculated the minimum and maximum clutch size used for the stretched b distribution (Morris and Doak 2002) as the mean clutch size ± (2 9 SD from the overall sample).
We used the number of females double-clutching (n) divided by the number of females telemetered during the reproductive period to estimate a mean d (d = probability of double-clutching). We calculated variance using the formula var(d) = d 9 (1-d)/n. We calculated the mean proportion and variance of double-clutching females for CWMA using the same method as that used for WR3.
We defined hatch rate as the proportion of hatchlings that successfully emerged from the nest (hatchlings divided by estimated clutch size). Successful eggs were easily distinguished from those that failed to hatch by appearance; white, papery eggs with a slit were assumed to be successful whereas dark brown, shriveled eggs with no evidence of an opening were considered failed. We estimated hatch rate for WR3 using field data from 2010 and 2011 and for CWMA using field data from 1998 to 2001. For each site, we calculated the hatch rate for each nest using an unweighted mean (because eggs in the same nest are not independent) over all nests (including first and second clutches). We used the variance of the mean hatch rate across all nests for the LSA.

Survival
We estimated annual survival rates of adult telemetered lizards that were tracked C10 days after marking during 2004-2011 on WR3 (n = 147 individuals) and during 1998-2005 on CWMA (n = 229 individuals). In a number of cases, we were unable to determine the fate of a lizard on WR3 because its transmitter signal disappeared, probably due to either transmitter failure or removal from the study area by a predator. Because of this ambiguity, we estimated survival rates in 2 ways that bracket the range of possibilities (defined as Categories; Endriss et al. 2007). Category 1 estimates assumed that lizards with undetermined fates were alive; these individuals were censored from the analysis following their disappearance. Category 2 estimates assumed that missing lizards were dead (Endriss et al. 2007). Lizards that died as a result of research activities were censored at the last date they were known to be alive. We used the Known-Fates model in Program MARK (White and Burnham 1999) to evaluate a priori hypotheses that sex, season (active or inactive), and study phase (for WR3) or year (for CWMA) affect survival, using the same time intervals and seasons as Wolf et al. (2013). Individuals tracked over multiple years were separated for the purposes of survival analyses, so experimental units for these analyses were individual-years. We used study phase, rather than year, for lizards on WR3 because of varying levels and types of disturbance to the site; these disturbances were described by Wolf et al. (2013).
Model-selection results, based on AIC c (Anderson et al. 2000), showed strong support for an effect of season in both populations. There was support for an effect of study phase on WR3 (Wolf et al. 2013), but no support for a year effect at CWMA. We constrained models to obtain overall mean and variance estimates for annual survival (accounting for season effects) at both sites to use for the LSA. For WR3 data, we averaged Category 1 and Category 2 estimates of the mean and variance of survival rates. There were no unknown fates in the CWMA dataset. We used the highest-ranked model including an effect of year (CWMA) or phase (WR3) to estimate variance among year/phase survival rates for input into the LSA. For CWMA, the estimated among-year variance combined with the low adult survival rate resulted in a distribution skewed towards 0 more than biologically likely based on survival analyses. We therefore reduced the variance from 0.028 to 0.020 to create a distribution with a mode closer to the mean and farther from 0. This adjustment did not change the coefficient of determination of k on adult survival by [0.02.
We calculated juvenile survival on WR3 for a given year, t, by dividing the number of juveniles recaptured as adults the following year, t ? 1, by the total number of juveniles marked in year t. We were unable to perform mark-recapture analyses because we had an inadequate sample of juveniles. However, Wolf et al. (2013) estimated a detection probability (p) of 0.4 (±0.07 SE) over 4-6 day intervals for adult P. cornutum on WR3. Extrapolating the detection rate for a 4-6-day interval over an entire 120-day active season by the equation 1-(1-p) x , where x = number of intervals, leads to a probability of detecting a given adult lizard (which may have been marked as a juvenile in the prior year) on WR3 of [0.99. Further, there was only one juvenile lizard recaptured in year t ? 2 and not in year t ? 1, and emigration appears to be minimal at WR3 (A. J. Wolf, unpublished data). We therefore feel that the detection issue is negligible when using the return rate to estimate juvenile survival. We averaged the survival rate across all years, and calculated an overall variance using the simple average of yearly variances. Data were missing from our geodatabase for 2003-2004 (recaptures in 2004 and 2005), and we used Endriss et al. (2007) to provide pooled data for the 2 years. Too few juveniles were recaptured to calculate juvenile survival for CWMA, and we therefore assumed it to be equal to the rate for WR3.
Detection and recapture rates of hatchlings were too low to provide an estimate of survival for either site. Moreover, very little is known about the causes of mortality in hatchling horned lizards. Both unimodal and bimodal distributions seem plausible, depending on whether main causes of death act on individuals (such as predators) or cohorts (such as climatic and meteorological effects like drought). Large variation in hatchling survival has been documented for other species of North American lizards in climatically variable ecoregions (Tinkle et al. 1993). Because of this uncertainty, we used a uniform probability distribution of hatchling survival for the LSA to include a wide range of possibilities. It is highly unlikely that hatchling survival is 1, so we used Euler's equation to estimate what hatchling survival rate would have to be for k = 1 given estimates of survival rates and fecundity for other stage classes (Hellgren et al. 2000;Endriss et al. 2007). We created a uniform distribution from 0 to twice the estimated rate of hatchling survival to obtain a plausible range of possible hatchling survival rates.

Model construction and analysis
Using the parameters and distributions from above, we populated a 3 9 3 Lefkovitch matrix (Crouse et al. 1987). This matrix was based on a life-history diagram that assumed 3 stage classes (hatchling, juvenile, and adult). In this model, we assumed that all hatchlings that survive grow to be juveniles, and likewise juveniles grow to be adults (i.e., the probability of surviving and not progressing to the next stage was 0 until adulthood). Adults surviving a year remained in the adult stage class the following year.
We used the PopTools extension (v. 3.2, G. Hood, Canberra, Australian Capital Territory, Australia) for Microsoft Excel to calculate k, deterministic elasticities, and perform Monte Carlo simulations. The Monte Carlo simulations drew 10,000 random values (Taylor et al. 2012) for each matrix element from their respective probability distributions, as well as calculated the resulting k for each set of randomly drawn matrix elements. We used PROC GLM in SAS 9.2 (SAS Institute Inc., Cary, NC, USA) to regress k against each set of randomly drawn vital rates, and calculated the coefficient of determination, R 2 , between each parameter and k (Wisdom and Mills 1997;Zar 1999). Several authors have suggested calculating covariance matrices for matrix elements in demographic populations, and these should ideally be incorporated into LSAs (Wisdom et al. 2000;Morris and Doak 2002;Johnson et al. 2010). However, this examination of covariance requires corresponding data over time (i.e., estimates of each vital rate for a number of years). Because many of the parameter estimates in this study were not available for multiple years (e.g., second clutch sizes at WR3 are only available from 2011), it was impossible to calculate these covariance matrices.
We conducted a perturbation analysis for hatchling survival because of uncertainty surrounding the probability distribution. We conducted LSAs for CWMA and WR3 with a variety of other minima and maxima for the uniform distribution (Biek et al. 2002;Taylor et al. 2012), encompassing a wide range of both means and sets of upper and lower bounds (which are analogous to variance in a b or normal distribution). Following preliminary regressions of k on hatchling survival, we selected 2 sets of bounds for the uniform distribution for hatchling survival in addition to the distribution based on the mean hatchling survival rate calculated from Euler's equation. We set one distribution such that both the minimum and maximum shifted slightly away from 0 and towards 1 (representing an increase in mean survival rate) and one with a smaller breadth (higher minimum and lower maximum, analogous to less variance). We ran LSAs and generated R 2 values from regressions of k on all vital rates for both populations using these 2 adjusted uniform distributions for hatchling survival rate. We report results for a subset of hatchling survival perturbations.
We also conducted a perturbation analysis for adult survival because the Category 1 estimates for adult survival on WR3 that we calculated in Program MARK were approximately double those calculated for Category 2 in each phase. To examine the effects of the assumptions involved in each Category of adult survival, we also conducted the LSA for WR3 using only Category 1 estimates, using only Category 2 estimates, using the mean survival rate of both Categories but the variance estimates from Category 1, and the mean of both survival rates but the variance estimate from Category 2. During these perturbations, we held all other vital rates constant. We used a lower bound of 0.05 and an upper bound of 0.65 for the hatchling survival distribution, thereby avoiding using hatchling distributions whose lower bound was 0.

Results
The vital-rate estimates for WR3 differed from those for CWMA in several aspects (Table 1). The adult survival rate at WR3 was more than double the rate at CWMA, whereas fecundity at WR3 was about 63 % of the corresponding estimate for CWMA. Clutch sizes at CWMA were larger than those at WR3, and its effect on fecundity was compounded by the greater proportion of females double-clutching at CWMA (Table 1). The deterministic elasticity of k to fecundity, hatchling survival, and juvenile survival was 0.26, and to adult survival was 0.23 at WR3. At CWMA, these respective elasticities were 0.31 and 0.08. Results of elasticity analysis of the deterministic matrix based on mean vital-rate estimates indicated that recruitment, including survival of the hatchling and juvenile stages and adult fecundity, had the greatest relative effect on k for both WR3 and CWMA. However, elasticity of k to adult survival was much closer to the elasticities of the other vital rates at WR3 than at CWMA, where it had less than one-third the effect on k compared to the other vital rates.
Results of perturbation analysis involving hatchling survival were similar for the 2 sites, and thus results from WR3 are not reported. Distributions of hatchling survival rates with a lower bound \0.01 led to the greatest effect of hatchling survival on k (Fig. 2). As the lower bound increased, the strength of the relationship between hatchling survival and k decreased incrementally (Fig. 2). However, if the minimum for the distribution was 0, the maximum of the distribution did not affect the R 2 values. For distributions with the same range between the minimum and maximum, those with upper bounds closer to 1 had lower R 2 values. Distributions with the same means but different ranges (e.g., 0.15-0.75 vs. 0.30-0.60) led to similar amounts of variation in k being explained. All subvital rates that composed fecundity had similarly low effects on k (Table 2). For WR3, nest survival, followed by hatch rate, consistently had the greatest effect on k of the vital rates that make up fecundity. Other subvital rates (proportion females reproductive, first clutch size, proportion females double clutching, and second clutch size) had very small effects on k. For CWMA, there were patterns of an increasing importance of subvital rates composing fecundity as minimum hatchling survival rates were increased, with hatch rate, and first and second clutch sizes being the most important parameters (Table 2).
Unlike hatchling survival, perturbations in the variance of adult survival had a greater effect on k than perturbations in the mean annual survival rate of adults for WR3 (Table 3). In all perturbations, hatchling survival still had the highest R 2 values, followed by juvenile survival. However, as the variance of adult survival increased, adult survival had greater effects on k to the point where it was approximately equal to fecundity. This effect was independent of the mean estimate used in the perturbation.
Overall, recruitment had a greater effect than adult survival on the growth of the population at WR3. Among the 3 versions of the LSA we ran for both areas, the version with hatchling survival with a lower bound of 0 explained the greatest variance in k (Table 2; Fig. 3). Under that scenario, the effects of hatchling survival on k as hatchling survival rates approached 0 dominated the other vital rates. As minimum hatchling survival rate increased, the effect of other vital rates increased and the effect of hatchling survival rates on k decreased accordingly, at least for CWMA. Juvenile survival had a greater effect on k than fecundity at WR3, whereas the effect of fecundity on k at CWMA was consistently greater than that of juvenile survival rate (Table 2; Fig. 3). The effect of adult survival on k in both areas was small (R 2 \ 0.1) under each scenario.

Discussion
Our prediction of a trade-off between adult survival and fecundity had some support, as the site with lower adult survival (CWMA) had higher fecundity and greater differences between the deterministic elasticities of k to adult survival and fecundity. In addition, adult survival had a greater relative importance at WR3 than CWMA, as predicted by the higher rate of annual survival at WR3. Despite the different importance of adult survival relative to fecundity at the two sites, however, deterministic elasticity analysis indicated that recruitment (fecundity, hatchling survival, and juvenile survival) were equally important at both sites, and outranked adult survival in importance. In addition, by adding stochasticity via the LSA, we did not find support for our prediction that we would observe population-specific contributions of individual vital rates to k. Although there were subtle differences between populations in the rank order of the importance of vital rates in affecting population growth, hatchling survival had the strongest relationship with population growth for these disparate populations of P. cornutum.
The higher adult survival at WR3 may be a result of a longer hibernation and a corresponding lower exposure to predators during activity (see Hellgren et al. 2010 andWolf et al. 2013 for seasonal survival rates). Conversely, a longer active season would allow lizards to eat and therefore grow more, leading to larger clutch sizes at more southern sites. Notably, hatchlings at both sites were the same size (R. T. Kazmaier and A. J. Wolf, unpublished data), so differences in hatchling survival based on initial size are unlikely. Studies of latitudinal variation in adult survival among lizards and snakes were consistent with our findings. Overwinter survival in side-blotched lizards (Uta stansburiana) increased with latitude across a 15°gradient (Wilson and Cooke 2004), and survival rates of adult eastern massasaugas (Sistrurus catenatus) increased along a southwestern to northeastern cline (Jones et al. 2012). Conversely, survival of black ratsnakes (Elaphe obsoleta) was higher at mid-latitudes (Sperry et al. 2010). From a theoretical perspective, our results were consistent with the hypothesis that among lizards, survival rates will be positively correlated to latitude, primarily via relationships among temperature, activity patterns, and associated predation risk (Adolph and Porter 1993), although other ecological factors may complicate this covariation (Sears 2005). In addition, survival and fecundity rates are also negatively correlated among populations occurring in different thermal environments (Adolph and Porter 1993). Recent work with Xenosaurus grandis, a lizard with higher survival among all stages and lower fecundity than Phrynosoma spp. (Zuniga-Vega et al. 2007), also documented the important influence of recruitment on k. For X. grandis, transition of hatchlings and juveniles to the next stage, as well as adult survival, had the greatest effects on k (Zuniga-Vega et al. 2007). Our results similarly highlighted the importance of recruitment, despite differences in life history between P. cornutum and X. grandis. Simulated P. cornutum population growth rates were most strongly related to recruitment in both study populations. In addition, although the P. cornutum populations differed in the rank orders of fecundity and juvenile survival, adult survival had the weakest effect on k.
A trade-off between juvenile survival and fecundity was observed in the demographic parameters of boreal toads (Bufo boreas; Biek et al. 2002), which exhibited lower juvenile survival than two frog species (Rana spp.), but had much higher fecundity. We did not have an independent estimate of juvenile survival for CWMA, but lower adult survival rates often correlate to lower juvenile survival rates (Pike et al. 2008). At WR3, juvenile survival had a stronger influence on k than fecundity. An estimate of juvenile survival at CWMA would reveal whether this holds true across populations of P. cornutum.
Hatchling survival was the top-ranked parameter in its effect on k at both sites, and this result was robust to perturbations of hatchling survival and adult survival (Tables 2, 3). Although the inclusion of 0 in hatchling survival distributions depressed the effects of all other parameters, a hatchling survival rate of 0 or a uniform distribution of rates among years are plausible assumptions. Highly variable survival rates of hatchlings have been documented in lizards (Heulin et al. 1997;Tinkle et al. 1993;Zuniga-Vega et al. 2007) and other reptiles (Hyslop et al. 2012;Perez-Heydrich et al. 2012). For example, hatchling survival of Scleroporus graciosus can vary annually from 0.12 to 0.59 (Tinkle et al. 1993). In addition, a meta-analysis of hatchling survival rates in the gopher tortoise (Gopherus polyphemus) concluded that a summary survival rate (±95 % CI) of 0.13 (0.04-0.34) should be used to parameterize future population models for that species (Perez-Heydrich et al. 2012). As they often hatch at the height of summer in hot, dry conditions, years of extreme weather could conceivably kill an entire cohort of hatchlings. It is of note to managers that years of no recruitment will negatively affect population growth in species with life histories similar to P. cornutum, regardless of adult survival and fecundity.
We acknowledge that the strength of the relationship between hatchling survival and k could be a result of our uncertainty in our estimate as much as it is a biological factor for population growth. Furthermore, we note that our model for CWMA was missing two of four vital rates (hatchling and juvenile survival). We stress the importance of accurately estimating these vital rates, about which very little is known for the order Reptilia (Pike et al. 2008). Potential techniques to monitor younger age classes in small (\100 g) herpetofauna include harmonic radar (Riley et al. 1996;Moseley et al. 2004) and genetic tracking. Continued efforts to develop these and other methods are critical for understanding reptile population dynamics.
Highly variable vital rates can have large effects on k (Biek et al. 2002). However, in our models, uniform distributions for hatchling survival rate that had the same mean but different upper and lower bounds (analogous to different variances for a b distribution) had similar R 2 values. For example, assuming a small minimal proportion (e.g., 0.05-0.15) of hatchlings survive did not change the rank-order results of the LSAs (Fig. 3). These results indicate that the tendency of a vital rate to take on extreme values stochastically, as opposed to its variance per se, may be causing the large contributions to variance in k. In other words, it is not only the variance that affects the importance of a vital rate, but also the likelihood of values very close to 1 or to 0. This finding bears noting in future LSAs.

Management implications
Our results have implications for conservation of threatened and endangered reptiles. First, they emphasize the importance of studying survival of non-reproductive stages. Due to the relative ease of capture and radiomonitoring adult Phrynosoma spp., the adult cohort is more commonly studied. This bias holds true for most vertebrates, including reptiles, with information on non-reproductive stages grossly lacking (Pike et al. 2008). Studies of adult survival, however, may be largely wasted if recruitment is driving population dynamics; it is therefore vital that managers identify cause-specific mortalities for juvenile and hatchling horned lizards, as well as other reptile species. However, it should be noted that subvital rates composing fecundity can only be quantified by studying reproductive females (which lead researchers to finding nests, an often difficult task with many reptile species).
Elasticity analysis and LSAs for our 2 disparate populations of P. cornutum indicated that augmenting nonreproductive stage-classes could increase population growth rates in P. cornutum and other species with similar life histories. Head-start programs, in which young animals are born and partially raised in captivity before release into the wild, have been attempted for a variety of reptiles and amphibians (Dodd and Siegel 1991;Ramo et al. 1992;Heppell and Crowder 1998;Spinks et al. 2003;Sprankle 2008). Generally, long-lived species (those with high adult survival) do not benefit substantially from head-start programs, as improving survival of young stages does not adequately compensate for adult mortality if adult survival is the key vital rate in these species (Crouse et al. 1987;Heppell and Crowder 1998;Enneson and Litzgus 2008). However, there are cases of successful head-start programs among species such as Chiricahua leopard frogs (Lithobates chiricahuensis; Sprankle 2008), and a variety of Iguana spp. (Escobar et al. 2010). Our demographic modeling demonstrates how head-start programs can be effective in species with low adult survival and high fecundity, such as P. cornutum.
The importance of fecundity and its components (nest survival, hatch rate, clutch size) in affecting k, although less than hatchling survival, are of interest to management. Minimizing human and habitat disturbances at a site may decrease the risk of depredation for nests (and lizards). Depredation of adults has increased in recent years at WR3 following human-induced habitat disturbance (Wolf et al. 2013), or adult survival in the absence of humans could be even higher at WR3 compared to CWMA. Clutch size and hatchling health upon emergence may be affected by female nutritional condition and body size, which may in turn be affected by habitat quality and food abundance (Ballinger 1983;Ford and Seigel 1989;James and Whitford 1994), factors that managers can conceivably modify. Unfortunately, hatch rate may be largely determined by climatic conditions. Attrition of hatchlings as a result of extreme climatic conditions (e.g., extended hot, dry weather) may be unavoidable for managers (with the exception of collecting and incubating eggs in captivity).
Whereas Johnson et al. (2010) found that differing vital rates among conspecific populations could lead to population-specific management strategies, our results suggest the opposite for P. cornutum. Although the absolute strength of the effect of each demographic parameter on k differed between horned lizard populations, management strategies focusing on hatchling survival would strongly benefit both populations. As a consequence, managers across the range of species with similar demographics to P. cornutum could adopt similar management priorities with respect to stage classes, despite intra-population differences in population vital rates. Biek et al. (2002) and Johnson et al. (2010) suggested that widely varying vital rates may drive k more than the vital rate ranked highest by deterministic elasticity analysis. Our results do not disagree with this conclusion, although whether variance itself or how frequently a vital rate stochastically approaches extremes (as a result of mean and variance) affects k more strongly needs further exploration.
In the wider arena of sensitivity analyses and conservation, quality demographic data (means and variances) for all stage classes must be obtained for wildlife populations of interest, as focus on life stages that are most easily monitored can risk missing the key factors that are driving population growth rates. Our study shows that drastic differences in vital rates between populations do not always require different management strategies, although the strength of our results is tempered by our inability to estimate survival rates of juveniles at CWMA and hatchling survival at both sites. Ultimately, our ability to fully understand population growth is linked to our ability to accurately estimate vital rates.