Local temporal trajectories explain population-level responses to climate change in saguaro ( Carnegiea gigantea )

. Population demography is typically assumed to be strongly in ﬂ uenced by climatic factors, particularly with succulent plants and cacti. The saguaro cactus ( Carnegiea gigantea ) is a long-lived columnar cactus of the Sonoran Desert that experiences episodic recruitment and mortality. Previous studies have attributed long-term changes in saguaro populations to climatic factors, including increased germination and establishment during wet periods and mortality and reduced establishment during droughts and extreme freezes. We used a 48-yr data set of marked individuals at the Desert Laboratory in Tucson, Arizona, to test the hypothesis that local, temporal population trajectories are mediated by topographic heterogeneity that interacts with ﬂ uctuating climatic conditions. We tested the in ﬂ uence of local slope and aspect vs. climatic variability on a population of saguaro using > 5800 marked individuals that have been measured since 1964. We examined the relationship between demography and climatic variables (drought, precipitation, and extreme temperatures) and found signi ﬁ cant differences in growth and survival among aspects and among census periods. Saguaro population growth was higher during wet and cool periods (e


INTRODUCTION
Spatio-temporal heterogeneity of biotic and abiotic factors plays an important role in modulating how climatic fluctuations affect population dynamics at local and regional scales (Lenoir et al. 2013).In plants, spatial heterogeneity buffers populations from fluctuations in environmental conditions that enhance or hinder plant establishment (Randin et al. 2009, Butterfield et al. 2010), and can promote long-term population persistence and species coexistence (Harrison and Taylor 1997, Chesson 2000, Sears and Chesson 2007).Factors producing spatial heterogeneity relevant for individual plant recruitment include physiographic factors such as landscape position, aspect, and slope among others (Turner 1990, Dobrowski 2011, Winkler et al. 2018).The effects of these factors can vary in space and time, interact with climatic variability, and influence species persistence differentially at various spatial scales or for a certain time (Graff and Aguilar 2017).For example, microhabitats that are suitable for the recruitment and persistence of a species oftentimes create a spatial signal in vital rates that are only detectable at certain spatial scales (Wisdom et al. 2000) and at the temporal scale of the species' demography (Dobrowski 2011).This is especially true in long-lived species with life histories that are difficult to measure and slow to respond to climate variability (Goldberg and Turner 1986, Butterfield et al. 2010, Kroiss and HilleRisLambers 2015).Thus, the uncertainties around which factors produce relevant spatial heterogeneities create difficulties when attempting to predict when and where plant species will respond to climate change.
Modeling populations of long-lived perennials is challenging given the time scale at which they respond to change (Morris et al. 2008) in addition to the spatial variation in the multitude of biotic and abiotic factors that influence population dynamics (Wright et al. 2014, Kroiss andHilleRisLambers 2015).The interaction of spatial and temporal heterogeneity adds further complication to understanding past and predicting future dynamics (Freestone and Inouye 2006, Conlisk et al. 2018, Winkler et al. 2019a, 2019b).
When key factors that benefit species persistence in otherwise unfavorable conditions are ignored, model predictions will likely underestimate species performance; likewise, models will overestimate species performance if they ignore factors restricting species persistence (Ashcroft 2010).For example, if precipitation and topography interact to produce spatial heterogeneity in water availability relevant for plant survival, predicting future species distributions based only on the independent effects of both factors will give wrong answers.As a result, the interactions of the full suite of climatic and topographic drivers of demographic rates need to be considered when modeling past and future dynamics.
Methods that aim to predict species responses to climate usually incorporate spatial heterogeneity as a modifier of response functions to climatic variables, but typically do not include interaction functions, though this is rapidly changing in recent years (Dobrowski 2011, Clark et al. 2014, Oldfather and Ackerly 2019).Environmentally induced variation in life-history traits such as survival rates, reproduction, or recruitment is typically modeled as random fluctuations of mean vital rates (Akcakaya et al. 2004), though this is also rapidly changing (Bellier et al. 2018, Yang et al. 2019).Additionally, spatial heterogeneity is usually treated as additive to climate by specifying these functions for different habitats.This approach largely ignores the interaction between climate and spatial heterogeneity.Recent studies have suggested that fine-scale topographic heterogeneity (at a scale of 0.01-1 km) in climate (topoclimate; sensu Ackerly et al. 2010) can significantly buffer populations from the local impacts of climate change (Randin et al. 2009, Willis and Bhagwat 2009, Ackerly et al. 2010).Thus, explicitly considering finescale topographic heterogeneity appears essential to predict species responses to future environmental change.Yet, a thorough approach that detects and explicitly incorporates correlations among life-history parameters and topographic variables is likely needed to accurately estimate changes in species' life-history transitions across a landscape.
The saguaro cactus (Carnegiea gigantea [Engelm.]Britt.and Rose) is an iconic columnar cactus in the Sonoran Desert with a range extending from southern Arizona south through Sonora and into northwestern Sinaloa, Mexico (Hastings 1961, Turner et al. 1995).Its ecological prominence in the upper Sonoran Desert and its importance as a food and habitat resource for many species underscore its importance as a keystone species (Steenbergh 1983, Turner et al. 1995).Numerous studies have addressed the ecology and growth rates of this species, mostly focusing on ageheight relationships to develop an expected population-level age structure, and then comparing that age structure to long-term climatic fluctuations (Pierson et al. 2013).Additionally, previous studies have shown that survival of small saguaro is reduced by freezing (Shreve and Hinckley 1937, Steenbergh and Lowe 1977, Nobel et al. 1986) and that individual growth rates and branch production (a proxy for reproductive potential) are favored under wet conditions during summer and winter (Drezner 2003, 2005, Pierson et al. 2013).Saguaro populations have episodic recruitment that coincides with multiyear periods of higher-than-average precipitation (Jordan and Nobel 1982, Turner 1990, Winkler et al. 2018).Some studies have shown that recruitment also relates to below-average maximum summer temperatures (Steenbergh andLowe 1977, Drezner andBalling 2002).This tight climatic control on recruitment and growth explains why the distribution of saguaro is limited by freezing temperatures in the northern part of its range and at higher elevations, and by reduced and unreliable summer precipitation in the western Sonoran Desert (Drezner 2005, Pierson et al. 2013).However, a lack of regional synchroneity in age structure suggests that regeneration might also be controlled by local biotic and abiotic heterogeneity (Pierson et al. 2013, Winkler et al. 2018), not just regional climatic fluctuations.
Additional studies on saguaro have found variable population-level responses to climatic fluctuations, which suggest that local conditions modulate population structure (Pierson et al. 2013).Recent work has also shown that local site conditions and landscape position (e.g., bajada, foothill, steep slope) can drive population recruitment in saguaro (Winkler et al. 2018).Within the northern Sonoran Desert, saguaro survival is greater on east-and south-facing slopes (Niering et al. 1963, Turner et al. 1995, Pierson and Turner 1998) because these aspects tend to experience less extreme summer and winter temperatures (Haase 1970), and less intense and shorter periods of freezing temperatures (Jordan and Nobel 1982).If local populations have unique trajectories in response to climatic fluctuations, predictions of regional or desert-wide changes in saguaro populations will likely be a challenge.However, it remains unknown whether local abiotic conditions and climate fluctuations independently influence population trajectories, or whether they interact with climatic fluctuations to produce unexpected regeneration patterns.Disentangling these interactions, or lack thereof, may be the important next step that synthesizes decades of saguaro research by revealing fine-scale controls on population demographics.Further, uncovering these patterns is critical to understanding past and current saguaro population dynamics, and to predicting saguaro and other long-lived species' responses to future climate change.
Here, we use a 48-yr-long, spatially explicit data set to assess the importance of local spatial heterogeneity on demographic responses to temporal climate variability in saguaro.We sought to answer the following questions: (1) "What do age structure and recruitment year of saguaro populations reveal about the sensitivity of the species to climate?" and (2) "How do growth, survival, and recruitment rates vary across habitats and time periods with distinctive climatic characteristics?"We expected that populations would recruit primarily during favorable conditions given the strong climatic relationships reported in previous studies.We also expected that saguaro populations would respond to topography at small scales but that these differences would be subtle and would not override the stronger climatic relationships.

Site description and measurements
We used long-term saguaro data from the Desert Laboratory on Tumamoc Hill (32°13 0 12.281″ N, 111°0 0 16.098″ W) in Tucson, Arizona (Rodr ıguez-Buritic a et al. 2013).The slopes on Tumamoc Hill range from gentle on the lower north side to steep along the south and east sides (Table 1).Saguaro at Tumamoc Hill has been individually monitored since 1964 when plots were established on north, south, east, and west slopes representative of the conditions of each aspect (similar sunlight exposures and substrates).Plots were established as a 250 m wide rectangle extending from the top to the bottom of the hill and are between 10.1 and 15.5 ha in area (Table 1).

Climate data
We calculated the annual number days with extreme temperatures (freezing or above 38°C) to identify periods that may explain changes in saguaro vital rates.We used linear regression to obtain slope estimates as a metric to quantify increases or decreases in the number of days with extreme events each year.In order to incorporate variation due to climate periods, we first categorized each census interval as either wet and cool, mild, or dry and hot using Palmer's Drought Severity Index (PDSI) for summer and cool-season precipitation separately.Palmer's Drought Severity Index combines precipitation and temperature influences on evapotranspiration into an index of aridity (lower values are more arid) that has been used in previous saguaro studies to explain regeneration patterns (Pierson and Turner 1998, Drezner and Balling 2002, Winkler et al. 2018).This classification was based on deviations from multiyear averages calculated from 1964 to 2012.We calculated PDSI in two ways: as an estimated PDSI using tree ring reconstructions from 1760 to 2004 and as an instrumental PDSI using reconstructions from 1895 to 2012 (see Appendix S1 for a description of calculations).We then used historical reconstructions of PDSI and coolseason precipitation to identify similar cool and wet, mild, warm, and dry periods between 1700 and 1964.Periods with climate that could not be classified within our system were classified as average (Appendix S1: Table S2).Each census period was classified differently based on its precipitation and PDSI values: 1964precipitation and PDSI values: -1970precipitation and PDSI values: (cool and wet), 1970precipitation and PDSI values: -1993precipitation and PDSI values: (mild), 1993precipitation and PDSI values: -2012precipitation and PDSI values: (warm and dry), and 1964precipitation and PDSI values: -2012 (average) (average).

Demographic vital rates and population structure
We summarized demographic vital rates for the entire population and each plot.First, we calculated growth rates (cmÁcm À1 Áyr À1 ) between each census for the entire population.We also calculated the range of growth rates for the entire sampling period .We then examined annual mortality (%) and the number of recruits for each plot as well as the entire population between each census period.Next, we tested for differences in growth rates across plots using ANOVA.Recruitment and mortality rates were analyzed across plots and census periods.
We then analyzed regeneration patterns of saguaro by comparing the observed and expected age structures for the populations in 1964 and 2012.Observed age structures were determined from age-height relationships based on field data (Turner 1990, Pierson and Turner 1998, Pierson et al. 2013).Expected age structures were derived as a constant vital rate prediction using age-specific mortality functions described in Appendix S2 (sensu Pierson and Turner 1998).Functions were scaled to the observed population density by multiplying the observed age structure and mortality to estimate the number of individuals expected in each height class during a given census period.This metric yields the population structure expected when recruitment compensates mortality and population structure is stable.Differences between these two age structures provide estimates of when cohorts were larger (positive residuals) or smaller (negative residuals) than required to maintain stable populations.Estimating age structures this way allowed us to account for three sources of variability: time periods with different climates, intraspecific variation in growth rates, and other unexplained sources of variation (i.e., error; Appendix S2).Notes: Plots were established in 1964 as 250 m wide rectangles that extend from the top to the bottom at each of the north, east, south, and west slope aspects of Tumamoc Hill.The north plot was divided into two subplots (lower and upper) to avoid disturbed areas around Desert Laboratory administrative buildings.
We then compared the observed age structures for each plot-time period combination with the expected age structure.First, we used plot-specific logistic-cubic-spline regressions to model the probability of mortality in a given year as a function of height (m) at the beginning of the census period (Appendix S2: Fig. S2).We averaged these predictions to produce 5-yr interval of survival (between 0 and 250 yr of age), which can be assumed to be the stable, expected age structure when mortality compensates recruitment.The final age structure obtained was scaled so the area under the curve was equal to the number of individuals in the population at each census (1964 and 2012) and within each plot in order to avoid overestimating age structures.
The difference between observed and expected age structures was used to quantify how a given cohort contributes to maintaining stable populations.Although observed age structures are related to both recruitment and mortality, we assume that residuals are a good approximation of past recruitment events since mortality can be described as a monotonic, decreasing function of age (or height) and random age-specific mortality events are likely rare (Drezner 2006a).Thus, positive residuals indicate more individuals than expected to compensate for age-specific mortality and we interpret them as recruitment peaks.We obtained normalized residuals first for the whole population on Tumamoc Hill and then separately for each plot.We did so using estimated age structures derived from lower-and upper-limit age estimates for the whole population.For each plot, we averaged lower-and upper-limit age estimates before obtaining normalized residuals.

Saguaro's relationship to climate
To explore the relationship between each cohort size and climate, we used Kendall's Tau Rank correlation coefficient (s) between normalized residuals and time series data for our selected climatic variables.Correlations were calculated for each time point from 20 yr prior to 20 yr after estimated germination dates.Climatic variables-estimated and instrumental PDSI, summer, and cool-season precipitation-were obtained from climatological reconstructions derived from regional dendrochronological data (Ni et al. 2002).Data for each variable were available for different periods, and data for all variables were only available after 1907.For this reason, correlations were estimated using all climatic variables and residuals post-1909, and then using only climatic variables with a longer time span and residuals since 1848 (i.e., the estimated year of the first cohort in all plots).Appendix S1 provides details on compilation and analyses of climatic variables for this analysis.Mortality gradually reduces the number of individuals remaining as cohorts age, thereby reducing our ability to detect significant correlations with climatic variables (Moorcroft et al. 2001).To reduce this source of bias, we dimmed climatic time series by applying a year-specific decay function that was the same mortality function estimated with data between 1964 and 2012.We calculated s with this decay-adjusted time series.All analyses were conducted in R version 3.0.2(R Core Team 2014).

Climate data
For the combined data set, the annual number of days with extreme temperatures saw a substantial decline beginning in the early 1900s (Appendix S1: Fig. S1).However, starting in the 1920s, the number of days with temperatures above 38°C increased by ~1 d every year (slope of the number of days vs. year = 0.52 d/yr, standard error [SE] = 0.05, P < 0.0001) throughout the period analyzed, with substantial changes occurring between 2002 and 2012 from 62.2 (SE = 6.4) to 95 d (SE = 12).After little change during the 20th century, the number of days with freezing temperatures increased by an average of 1 d every 10 yr (slope of 0.1 d/yr, SE = 0.18, P = 0.045) with a sharp increase from 16.1 (SE = 8.8) to 48.3 d (SE = 24.3) between 2002 and 2012 which drove the relationship.No other climate variable had a significant long-term trend.

Demographic vital rates and population structure
The demographic rates for saguaro populations on Tumamoc Hill showed significant variation among census periods and among plots (Table 2).Both annual mortality and growth rates declined from 1964 to 2012, while annual recruitment was lowest for the mild middle period of 1970-1993 (Fig. 1).The south plot had higher mortality between 1964 and 1993 (which is classified as cool-wet and mild in climate) than the other three plots which saw similar amounts of mortality.The difference was mostly due to a slightly higher mortality rate of intermediateheight individuals in the south plot (Appendix S2).
Annual mortality of all plots was highest from 1964 to 1970, a period dominated by cool and wet climatic conditions (Appendix S1: Table S2).The lowest mortality occurred during the warm and dry period of 1993-2012.Across all plots, mortality was higher for very small and very Notes: Mortality percentage is calculated within each period/plot combination without considering new recruits.Mortality is divided by the census period to give the annual mortality reported.New recruits include new plants at the end of the census period and previously un-detected plants.
† The reproductive potential at the end of each census estimated as the percentage of plants presumed tall enough to be reproductive (≥2 m) out of the total estimated population size; the number of potentially reproductive meristems (number of individuals plus number of branches) is given in parenthesis.
‡ Total estimated population size includes the detected individuals at the end of each census periods plus those plants that germinated before but were not detected.large individuals.Two exceptions to this trend were higher rates of mortality for small-to intermediate-sized individuals in the south plot in the mild period between 1970 and 1993, and for intermediate-and large-sized individuals in the north plot in the dry and warm period between 1993 and 2012 (Appendix S2).
Recruitment rates estimated with normalized residuals from 1964-and 2012-observed age structures showed a long-term pattern that was congruent between lower-and upper-limit estimates (Fig. 2, top panels).Residuals from data generated using 1964 age structures showed a larger than expected recruitment during the first third of the 20th century, with either a wet or an average climate, and lower than expect recruitment before or during the early years of the mid-century drought of 1946-1950.A large recruitment event occurred during the second half of the century, but it was inconsistent between lower-and upper-limit estimates.Estimates using data generated using 2012 age structures showed more consistency for this second recruitment peak around 1950.Negative Fig. 2. The normalized residuals between estimated and expected age structures for the entire and plot-specific saguaro populations at Tumamoc Hill.Residuals for (a) 1964 and (b) 2012 (right); lines represent residuals using different age estimation limits (solid black = upper limit estimations, solid gray = lower limit estimations); time on the x-axis is divided into periods of contrasting climatic conditions (vertical dash-dot-dash lines) represented by a, average; w, wet; m, mild; and d, dry (see Appendix S1 for details regarding these estimations); black dots mark when censuses were conducted (1964, 1970, and 2012).(c-f) residuals for each plot, averaged between upper and lower residuals for 1964 (gray solid lines) and 2012 (black solid lines), with time range divided into periods of contrasting climatic conditions.

❖
8 August 2019 ❖ Volume 10(8) ❖ Article e02844 recruitment residuals were consistently detected from 1978 onward.Differences between lowerand upper-limit estimates in how strong peaks of recruitment were and when they occurred reflected the uncertainty around age estimations, while differences between data sets (1964 vs. 2012) primarily reflected mortality, suggesting that mortality events between 1964 and 2012 did not change the overall population structure.
The effect of aspect represented by plot-level differences mostly corresponded to discrepancies with when recruitment events occurred or how large they were.East and south plots tended to behave alike, with large recruitment events occurring in the first half of the 20th century, especially during the wet period of 1910-1922, and at the beginning of the drought between 1946 and 1950.After 1964, large recruitment events occurred during the dry period of 1971-1977 (Fig. 2e-f, Fig. S1b).In contrast to these plots, a recruitment pulse ca. the 1950s was small or missing in the north and west plots, and large recruitment events were more frequent but not as large, with peaks between 1861 and 1945 for the north plot and between 1826 and 1945 for the west plot (Fig. 2c-d).We found a unique pattern in the north plot: The population along this aspect changed from initially being dominated by relatively old individuals to being dominated by individuals <30 yr old in 2012 due to a huge late 20th century/early 21st century spike in recruitment.This change likely explains the inconsistent presence of the 1850-1922 recruitment peaks (positive residuals) between data generated using 1964 and 2012 age structures (Fig. 2c).In the last 47 yr, only the north plot experienced higher than expected recruitment events, and it was also the only plot with a large change in population structure (Fig. 2c).

Saguaro's relationship to climate
Contrary to our expectations, we did not find an indisputable pattern of positive correlations between large recruitment events and favorable growing conditions for saguaro (i.e., low drought, high precipitation, few days with temperatures below zero or days >38°C; Fig. 3).Residuals based on 1964 age structures revealed large recruitment events correlated with periods of low drought severity (high PDSI values) and high precipitation for summer and winter for at least some lag values.Small recruitment events (negative residuals) increased with number of days above 38°C (in Fig. 3, s > 0.1 or <À0.1 corresponds to P < 0.05) with this earlier data set.For residuals estimated with 2012 age structures, which adds dynamics for the last 48 yr to the data set, large recruitment events correlated with poor growing conditions before and after estimated germination times (Fig. 3a, b).These patterns hold when using residuals and climatic time series from 1909 or from 1848 to 2012 (Fig. 3a-c; Appendix S3: Fig. S1).The only significant difference between results for these two time series is the stronger correlation of recruitment with high cool-season rain found with the 1848-2012 time series.This suggests that larger recruitment events may have occurred more frequently during periods of high cool-season precipitation before 1909 than during the last century (s > 0.1 at most lags for south and east plots; Fig. 3c; Appendix S3: S1).Overall, this suggests that saguaro's response to climate in the second half of the 19th century was different from that after the start of the 20th century.
Plot-specific correlations show few commonalities (Fig. 4) which suggest a strong effect of aspect on population response to climate.Only two trends at the population level consistently downscaled to the plot level.Large recruitment events positively correlated with low drought severity and high cool-season precipitation before 1964 (Fig. 4a-c), while correlations became weaker when residuals up to 2012 were used.In addition, large recruitment events were mostly correlated with unfavorable climatic conditions when these newer (up to 2012) data were analyzed.Aside from these two trends, aspect plays a role in which and how climatic variables are associated with large recruitment events.
There were three distinct aspect-specific patterns.First, the north plot showed different trends from the other plots for most correlations.Namely, large recruitment events in the north plot were consistently correlated with favorable PDSI and winter precipitation.However, correlations of recruitment with summer precipitation and number of extreme warm days were negative using the data sets generated with 1964 age structures but positive using the data set generated with 2012 age structures.This suggests favorable summer precipitation may have compensated for the more extreme high temperatures in the last 48 yr in allowing for high recruitment in the north plot.Overall, the north plot was the only one that showed a negative correlation between large recruitment events and number of freezing days for lag times near zero (Fig. 4e), suggesting a heightened vulnerability to cold weather.
The second pattern is that the west plot had similar trends, but different correlation values compared to south and east plots.This was consistent using residuals from age structures from both 1964 and 2012 and for all variables except number of freezing days.The correlations with PDSI, precipitation, and warm days before 1964 were stronger for the west plot (Fig. 4a, b and e, f), while correlations with PDSI and precipitation using the 2012 ages structure were insignificant for the west plot (Fig. 4c-d).This suggests that the west plot was the least sensitive to variation in climatic conditions during the last 48 yr.
The last pattern was that the east and south plots had similar recruitment trends that showed similar relationships with climatic variables (Figs. 2, 4).Before 1964, the sign of the correlations varied from positive to negative at different lags for all variables, indicating that results depend critically on lag assumptions.Using the 2012 residuals, recruitment was mostly correlated with years of drought (low PDSI), low precipitation, and more freezing days; conditions usually considered to be unfavorable for saguaro recruitment.The only favorable correlation using the 2012 data was the association of recruitment with years with fewer extremely warm days.

DISCUSSION
Our study of saguaro demonstrates how vital rates at local scales track seasonal precipitation, drought severity, and temperature to create subpopulations with varied demographics across our study site at the Desert Laboratory in Tucson, Arizona.We confirmed that overall recruitment patterns in saguaro are in accord with previously proposed controls of population-level dynamics, but also found unexpected responses under contemporary climate.First, like previous studies, we found that prior to 1964 recruitment was more likely to occur during favorable climatic conditions that include low drought stress (Jordan and Nobel 1982, Turner 1990, Winkler et al. 2018), above-average summer precipitation, and below-average maximum summer temperatures (Steenbergh andLowe 1977, Drezner andBalling 2002).However, after 1964, recruitment was more likely under unfavorable conditions (i.e., high drought severity, low winter precipitation).This is especially surprising in the north plot where a large recruitment spike occurred under high drought severity and in years with many extremely warm days.Our study demonstrates that, while recruitment is episodic, in the future, complex dynamics should be expected with a lower sensitivity of recruitment to the previously important limiting factors such as below-freezing temperatures (Steenbergh and Lowe 1977).
Population changes under fluctuating climatic conditions often result from a complex interaction of life history and climatic stresses in which sensitivity to changes strongly affects long-term trajectories (Oostermeijer et al. 1994, Conver et al. 2017, Winkler et al. 2019b).As a result, we expect highly sensitive populations to respond more strongly and positively to periods that favor recruitment and survival but also strongly and negatively to unfavorable periods.We found a clear signal of the influence of aspect in mediating how sensitive saguaro populations are to climatic fluctuations.Our study is the first to show saguaro recruitment occurring nearly independently of climatic conditions (west plot, last 48 yr), concurrently with unfavorable conditions including intense drought and low precipitation (east and south plots), and under a combination of favorable-high summer precipitation-and unfavorable-high frequency of warm daysconditions (north plot).These clear response differences between our four plots demonstrate the importance of local topography in mediating long-term population responses.
We suggest that the interaction between topography and climate operates through multiple processes.First, climatic variables may compensate for one another on different aspects to favor recruitment pulses.This is exemplified by contemporary dynamics over the last 48 yr in three plots where recruitment was associated with unfavorable values of at least one climate attribute.Second, aspect also appears to mediate the influence of population density in buffering responses to climatic fluctuations.This is seen in our most densely populated plot (south) which had opposite responses to climatic attributes.Our south plot had the highest number of recruits but also the highest mortality rates.There was also a constant saturation of sites with potential new recruits despite fluctuating climate attributes and high mortality rates.Factors other than density might be particularly important along west and east plots given that we observed intermediate to high densities there but also found significant recruitment correlations with unfavorable conditions.
Density differences in saguaro across slopes have been explained by a lower probability of freezing-related deaths on the east and south slopes given that direct sun exposure comes sooner during sunrise, shortening freezing periods (Pierson and Turner 1998).We did not find strong evidence of this freeze protection in our analyses but instead found that only our north plot showed reduced recruitment in response to the number of extremely cold days (Fig. 4e, north plot at lag 2-8 using the 2012 correlation estimates).Although freeze-related mortality might be important for culling saguaro individuals (Nobel 1980, 1982, Drezner 2006a), over the long-term and at the temporal scale of our study, its negative effects on recruitment are likely canceled out by the positive effects of other factors creating recruitment peaks during favorable growing conditions.In addition to the effect of freezing, density differences among plots might result from differences in summer conditions.Drought conditions are likely reduced on the steep, rocky slopes of the south-and east-facing plots (Gitlin et al. 2006, Winkler et al. 2018), and under nurse plants such as palo verde (Parkinsonia microphylla; Steenbergh andLowe 1976, Drezner 2006c).However, with increasingly warmer winter temperatures and summer precipitation that have occurred in recent decades, seedlings may not benefit as much from the positive effect of rocks and nurse plants in areas previously at risk of freezing or drought-related mortality.
At Tumamoc Hill, recent recruitment has been insufficient to compensate for mortality.Saguaro populations are likely to continue to decline if this trend continues, confirming the predictions of a previous study that used the 1964-1993 data (Pierson and Turner 1998).This trend is not a result of low seed production given that reproductive potential of this population had not diminished by 1993 based on the number of flowering-sized plants (Pierson and Turner 1998) or by 2012 (Table 2).In the first decade of the 21st century, there has been a large increase in the number of days with extreme temperatures (both cold and hot; Appendix S1), and it is uncertain whether these climatic trends will continue.If they do, as is predicted by climatic modeling (Seager et al. 2007, Seager andVecchi 2010), our analyses suggest that saguaro populations on Tumamoc Hill will continue to decline.
Higher than expected recruitment on the north slopes does not ensure the persistence of saguaro on this plot, which had the lowest density of the four plots.The north plot also had the highest mortality of adults between 1964 and 2012, which shifted the age structure toward a juvenile-dominated population.This suggests a potential boom-bust population dynamic on the north slope that jeopardizes long-term population stability because low density increases vulnerability to unfavorable growing conditions (Wilmers et al. 2007).In addition, the south, east, and west plots have had insufficient recruitment to compensate for mortality during the period of measurement.Although their population sizes might be sufficient to buffer them from unfavorable growing conditions (Conver et al. 2017), the 21st-century trend toward more severe drought and extreme temperature events (Cook et al. 2004) might jeopardize the persistence of desert species even on these more favorable slopes.
Our results show that recruitment peaks are not perfectly synchronized with periods of favorable climatic conditions.We explicitly included sources of variation that were not included in previous studies (Turner 1990, Parker 1993, Pierson and Turner 1998, Drezner 2006b, Drezner and Balling 2008, Pierson et al. 2013).These included individual variation and a variable response of growth rates to changing climate conditions.When these sources of variation are included in our models, our age estimates had an average error of AE6 yr, which reduces our ability to detect population-level responses to episodic, annual events like high summer precipitation or extreme freezing temperatures.Instead, this level of resolution is best suited to assessing the importance of multiyear recruitment peaks over many decades.Moreover, our preponderance of low-correlation coefficients with climatic variables reinforces the conclusions of Pierson and Turner (1998) and Winkler et al. (2018) who found that other factors, such as microhabitat differences, are also important drivers of saguaro recruitment.Herbivory, competition with exotic grasses (red brome, Bromus madritensis ssp.rubens; and buffelgrass, Cenchrus ciliare), and nurse-plant effects are some of the biotic factors we did not analyze that are likely important for explaining additional plot-level and even microscale variability (Pierson andTurner 1998, Snyder andChesson 2004).Total cover by C. ciliare, estimated from high-resolution aerial photography, showed an increase of 170% between 1994 and 2000 and an additional 69% in the next seven years (J.L. Betancourt and T. Bean, unpublished data; Olsson et al. 2012), especially near the top of Tumamoc Hill on the south, east, and west slopes.Competition with exotic grasses has also been shown to influence saguaro range expansion potential (Springer et al. 2015, Swann et al. 2018).No other factors are known to correspond to the current reduced recruitment, although the role of increased mortality rates of P. microphylla (Turner et al. 1995) and other large perennials that serve as saguaro nurse plants (Breshears et al. 2005) are in need of further investigation.
Although we do not know whether changes in the saguaro population at Tumamoc Hill are typical throughout the species' distribution, this long-term study is one of the firsts to investigate decadal-scale demographic patterns in relation to local topographic differences.Our results suggest that predictions about a demographic response to climate require an explicit consideration of local site conditions, which are only detectable using long-term monitoring (Winkler et al. 2018(Winkler et al. , 2019)).Previous work that documented high saguaro recruitment on east-and south-facing slopes (Pierson and Turner 1998) would have been insufficient to detect the current population-level shifts we observed in response to climate.For example, north-facing slopes provide suitable recruitment niches for saguaro and the population there has responded with recent, above-average recruitment under warmer conditions in the early 21st century.Considering the current predictions of a higher frequency of dry years and warmer temperatures for the Sonoran Desert (Seager et al. 2007, Notaro et al. 2010, 2012), it remains to be seen to what extent the predicted northward shift in saguaro distribution (Notaro et al. 2012) also comes with topographic-driven changes in local populations toward northerly slopes.
Overall, our results have broad implications for managing the ecological consequences of climate change.Ecologists must integrate habitat distribution and population modeling to understand potential changes in species distribution ranges (Akcakaya et al. 2004, Araujo and Guisan 2006, Keith et al. 2008, Brook et al. 2009).This integration occurs by using bioclimatic models to predict future suitability of a given area and then using metapopulation models to apply this modified scenario.This also implies translating habitat suitability into changes in carrying capacity at every location.Current approaches only consider climate-driven changes in geographical range area or quality and quantity of suitable habitats (Araujo and Guisan 2006).A hypothetical habitat-distribution model for saguaro on Tumamoc Hill, based only on presence, presence/absence, or abundance data before 1993, would most likely have failed to predict the shift in recruitment patterns we observed, even when complemented with an understanding of physiological constraints.These shifts were only detected by analyzing expected-vs.-observedpopulation-age structures estimated from long-term monitoring.We agree with Conlisk et al. (2013), who suggest bioclimatic-envelope predictions could be improved by the inclusion of long-term demographic data given the sensitivity of the results to population model parameters.We also agree with Ibanez et al. (2013), who discuss a need for a more mechanistic understanding of specieslevel responses to climate.Given that long-term demographic data are scarce, our immediate challenges should be to define the minimum attributes a data set should include to precisely predict demographic response to climate change.

Fig. 1 .
Fig. 1.Vital rates at different census intervals for our saguaro population at Tumamoc Hill.(a) Annual growth rate in height (cm/yr) relative to the height at the beginning of each period for each of the three census periods and for the entire monitoring interval (1964-2012).Letters placed each boxplot indicate significantly different groups based on Tukey's multiple comparisons with a = 0.05.(b) The annual mortality percentage for the whole population and for each of the four plots (east, north, south, and west).(c) The annual number of recruits in each plot and for the whole population using the same line styles as in (b).

Fig. 3 .
Fig. 3. Correlations between climatic and normalized residual time series for saguaro at Tumamoc Hill using combined plot data.Each panel (a-f) depicts correlations between normalized residuals and one climatic variable at different lags (climatic variable values from 20 yr before to 20 yr after the year for the estimated residual).The legend in the top left panel (a) applies for all panels; black lines indicate estimated residuals based on 1964 age structures, and gray lines indicate estimated residuals based on 2012 age structures either using the whole range of residuals from 1848 to 2012 (dashed lines) or only residuals from 1909 to 2012 (solid lines).The solid horizontal lines represent the a = 0.05 significance cutoff.

Fig. 4 .
Fig. 4. Plot-specific correlations between climatic time series and normalized residuals after 1909.Each panel (a-f) depicts correlations between normalized residuals and one climatic variable at different lags (climatic variable values from 20 yr before to 20 yr after the year for the estimated residuals).The legend in the top left panel (a) applies for all panels; line color indicates which age structure was used during the correlation estimations (black = 1964 and gray = 2012).All correlations were calculated using residuals from 1909 onward (see Appendix S3: Fig. S1 for correlations with residuals starting from 1848).The solid horizontal lines represent the a = 0.05 significance cutoff.

Table 1 .
Topographic characteristics of the four saguaro (Carnegiea gigantea) permanent plots at Tumamoc Hill, Arizona, USA.

Table 2 .
Individual counts per census period at each of the four saguaro plots at Tumamoc Hill, Arizona, USA.