Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Changing Climate Drives Divergent and Nonlinear Shifts in Flowering Phenology across Elevations

Changing Climate Drives Divergent and Nonlinear Shifts in Flowering Phenology across Elevations Article Changing Climate Drives Divergent and Nonlinear Shifts in Flowering Phenology across Elevations Highlights Authors d Flowering time is diverging among communities across an Nicole E. Rafferty, Jeffrey M. Diez, elevational gradient C. David Bertelsen d Divergence reflects nonlinear shifts in flowering phenology Correspondence over three decades nicole.rafferty@ucr.edu d Climatic variables have also changed differently across In Brief space and over time Using 33 years of data, Rafferty et al. test d Changing climate is driving phenological reshuffling across whether flowering times have diverged local spatial gradients for neighboring communities comprising 590 species across an elevational gradient. Divergent and nonlinear shifts in flowering time are related to elevation- specific changes in temperature and precipitation and will likely alter species interactions and gene flow. Rafferty et al., 2020, Current Biology 30, 432–441 February 3, 2020 ª 2019 The Author(s). Published by Elsevier Ltd. https://doi.org/10.1016/j.cub.2019.11.071 Current Biology Article Changing Climate Drives Divergent and Nonlinear Shifts in Flowering Phenology across Elevations 1,2,6, 3 4,5 Nicole E. Rafferty, * Jeffrey M. Diez, and C. David Bertelsen Department of Evolution, Ecology, and Organismal Biology, University of California, Riverside, 900 University Avenue, Riverside, CA 92521, USA Rocky Mountain Biological Laboratory, PO Box 519, Crested Butte, CO 81224, USA Department of Botany and Plant Sciences, University of California, Riverside, 900 University Avenue, Riverside, CA 92521, USA School of Natural Resources and the Environment, University of Arizona, 1955 E. Sixth Street, Tucson, AZ 85721, USA Herbarium, University of Arizona, PO Box 210036, Tucson, AZ 85721, USA Lead Contact *Correspondence: nicole.rafferty@ucr.edu https://doi.org/10.1016/j.cub.2019.11.071 SUMMARY climate is driving phenological reshuffling across local spatial gradients. Climate change is known to affect regional weather patterns and phenology; however, we lack under- INTRODUCTION standing of how climate drives phenological change across local spatial gradients. This spatial variation is Flowering phenology is a key biological indicator of current climate change [1–5]. However, the magnitude and direction of critical for determining whether subpopulations changes in flowering times appear to vary significantly both and metacommunities are changing in unison or among species and among communities, in part because of vari- diverging in phenology. Divergent responses could ation in how much abiotic factors such as temperature and pre- reduce synchrony both within species (disrupting cipitation have changed and in how sensitive species are to gene flow among subpopulations) and among spe- those changes [6]. Variation in responses has also been linked cies (disrupting interspecific interactions in commu- to biotic factors. For example, in some communities, flowering nities). We also lack understanding of phenological phenology responses are related to plant traits, such as flower- change in environments where life history events ing season [5], whether species have an annual or perennial life are frequently aseasonal, such as the tropical, arid, cycle [2,5], and whether species are wind or animal pollinated and semi-arid ecosystems that cover vast areas. [5]. In general, species that flower early in the season or are an- Using a 33-year-long dataset spanning a 1,267-m nuals show the greatest advances in phenology, whereas evi- dence regarding pollination mode is mixed [2,5]. In some cases, semi-arid elevational gradient in the southwestern closely related species exhibit similar responses; thus, phyloge- United States, we test whether flowering phenology netic relationships can also be important predictors of shifts in diverged among subpopulations within species flowering time [7,8]. and among five communities comprising 590 spe- Altered flowering phenology in response to climate can affect cies. Applying circular statistics to test for changes population dynamics and demography via several avenues [9]. in year-round flowering, we show flowering has For example, plant species that track climate by advancing become earlier for all communities except at the phenologically have higher metrics of performance such as highest elevations. However, flowering times shifted individual growth [10]. Delayed flowering has been linked to at different rates across elevations likely because compression of the flowering period and lower fruit and seed of elevation-specific changes in temperature and set [11]. The timing of flowering in relation to snowmelt and precipitation, indicating diverging phenologies of damaging frost events can have large effects on floral abun- dance in natural populations of wildflowers [12], although conse- neighboring communities. Subpopulations of indi- quent floral abortion might not translate into reduced population vidual species also diverged at mid-elevation but viability [13]. Nonetheless, phenological responsiveness has converged in phenology at high elevation. These been shown to be an important target of natural selection [8]. changes in flowering phenology among communities Increased variability in flowering phenology among species and subpopulations are undetectable when data are over time has been detected in some temperate datasets [14] pooled across the gradient. Furthermore, we show and is also likely to have demographic consequences. For that nonlinear changes in flowering times over the example, greater variation in species’ flowering times could alter 33-year record are obscured by traditional calcula- temporal overlap among different species, potentially affecting tions of long-term trends. These findings reveal pollen transfer and seed set. greater spatiotemporal complexity in phenological While communities undergo temporal reshuffling, shifts responses than previously recognized and indicate in flowering phenology can alter interspecific interactions, as 432 Current Biology 30, 432–441, February 3, 2020 ª 2019 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). well. For example, earlier flowering can cause asynchrony be- variation in temperature [31], although key abiotic factors vary tween plants and pollinators, particularly for species that flower in different ways along elevational and latitudinal gradients in early spring and are visited by the earliest emerging pollina- [32]. Our primary goal is to determine whether flowering phenol- tors, such as queen bumble bees [15]. Experimental manipula- ogies of communities and subpopulations are shifting in the tions of phenology have demonstrated that flowering time can same ways across elevations or are responding differently over affect the frequency of both plant-herbivore [16] and plant-polli- space. Differential shifts at the community or subpopulation nator interactions [17,18] and how effectively plants are polli- levels would signal changes in temporal overlap or synchrony nated ([19], but also see [20]). Over longer timescales, phenolog- across elevations. A second goal is to determine whether chang- ical mismatches between plants and pollinators can account for ing temperature and precipitation patterns are responsible for the loss of some interactions and have likely contributed to local any community-specific changes in phenology. Since the late species extinctions [21], although in other cases, synchrony be- 1990s, the study area has been impacted by a long-term warm tween plants and pollinators is maintained under climate change drought that is characterized by not only precipitation deficits [22]. Indeed, analyses of pollination success over very long time- but also rising temperatures, increasing precipitation intensity, scales could be necessary to determine trends [23]. Competitive and increasing precipitation variability [33–36]. This drought and facilitative interactions between plant species are also likely and longer-term climatic changes could be affecting commu- to be reshaped by phenological shifts [24]. nities and subpopulations along this gradient differently, in part Overall, research to date has provided some understanding of because of species- and population-specific responses and what traits predict species-level flowering phenology shifts and spatial environmental variation that could modulate the changing how those shifts influence performance, demography, and spe- abiotic conditions. Thus, we seek to determine whether eleva- cies interactions. However, how long-term shifts in flowering tional differences in climate variables and their changes over phenology vary among species-rich communities and subpopu- time lead to temporally nonlinear shifts in community-level lations across semi-arid environmental gradients has not been phenology that could result in spatially divergent phenological previously investigated. Prior work on change in flowering changes. phenology across spatial gradients has instead addressed The dataset we analyze is unique; in addition to capturing a topics such as how shifts in flowering time vary among species steep elevational gradient comprising five communities (defined [3], across a mosaic of moisture habitats at a single elevation by elevation bands) with six intergrading vegetative associa- [25], or over only a few years [26]. Ultimately, divergence in flow- tions, it includes 590 vascular plant taxa ranging from annuals ering phenology among subpopulations (i.e., spatially structured to trees and was collected in a semi-arid ecosystem where pre- subsets of larger populations) is expected to result in decreased cipitation and temperature are key triggers of flowering for many pollen flow and greater reproductive isolation because conspe- species at low and high elevation, respectively [3,4,37]. There is cifics overlap less in flowering time across space [27,28]. little overlap in the flowering assemblages of the five commu- Equally, convergence in phenology via increased overlap in flow- nities except in winter [37]. Fifty percent or more of species in ering time across space could result in increased pollen flow and each community are highly opportunistic, wherein flowering is outbreeding depression if subpopulations are locally adapted triggered by antecedent climatic conditions [38–40]. Because [29]. At the metacommunity level, divergent phenological re- these species have one to several flowering periods of varying sponses among adjacent communities will likely affect trophic durations, and flowering late in one calendar year could be the and non-trophic interactions, leading to altered community early part of a flowering season continuing into the next calendar structure and ecosystem processes [30]. For example, if flower- year, we apply circular statistical methods to detect spatiotem- ing phenology in a montane plant community shifts unevenly poral patterns. across elevations, this will alter the timing of resource availability We address the following questions: (1) at the metacommunity for species in other trophic levels, such as pollinators, and will level, has flowering phenology diverged or converged across el- likely have downstream effects on seed and fruit production, evations; (2) can variation in temperature and precipitation affecting frugivores, plant recruitment, and competitive interac- across elevations explain shifts in flowering phenology; and, (3) tions. Despite the many possible ways in which climate- within species, has flowering become more, or less, synchronous change-driven shifts in flowering phenology could affect sub- for subpopulations across space? Our approach reveals spatially populations and metacommunities, investigation of these divergent, temporally nonlinear changes in flowering phenology topics is limited by a lack of long-term phenological data across among communities and subpopulations that are obscured by the spatial scales important for maintaining gene flow and the pooling of phenological data across space and the simple species interactions. Without such data, it is impossible to deter- linear calculation of long-term trends. These results highlight mine whether species are responding at finer spatial scales that the utility of circular statistics for detecting patterns in phenology could alter the temporal overlap among subpopulations and and point to a greater complexity of responses than has previ- communities. ously been recognized, suggesting that climate change will Here, we analyze a 33-year record of flowering phenology that reshuffle communities in multiple dimensions. spans an elevational gradient of 1,267 m in the southwestern United States that encompasses desert scrub at low elevation, RESULTS several semi-arid associations (riparian scrub, scrub grassland, oak woodland, and oak-pine woodland), and pine forest at Change in Flowering Phenology from 1984–2016 high elevation. This elevational gradient captures a range of The relationship between flowering time and year varied environmental conditions that could be indicative of latitudinal by elevation band, indicating elevation-specific phenological Current Biology 30, 432–441, February 3, 2020 433 Figure 1. Change in Flowering Phenology of Each Community (A–E) Circular plots show the dates of flowering, represented by colored lines, for two time periods for all species within each community ascending the transect. Within each circular plot, the lengths of the colored lines indicate the number of phenological observations on each date. The lighter color represents data from 1984–1993, and the darker color represents data from 2007–2016; numbers give the day of year on which each month begins. The arrows show the mean di- rection and the mean resultant length (a metric of concentration) for a given time period; shorter arrows indicate more dispersed flowering times. (A) 945–1,079 m (community 1); (B) 1,079–1,372 m (community 2); (C) 1,372–1,671 m (community 3); (D) 1,671–1,939 m (community 4); (E) 1,939–2,212 m (community 5). (F) The summary plot shows the change in flowering time for each community. See also Figure S1; Data S1A and S1B. responses of communities (Data S1A). We therefore report re- predictor; Data S2B). No significant change in flowering time sults per community, numbering communities from 1–5 and was detected for the highest-elevation community (5; Data S2A). ordering from the lowest (1) to highest (5) elevation bands. When examined with linear circular models, from 1984–2016, Change in Flowering Phenology from 1984–1993 versus communities 1–4 shifted significantly earlier at a rate of 2.5, 2007–2016 1.6, 0.44, and 0.36 days per year, respectively, and shifts An examination of changes in flowering time during the first and became smaller with increasing elevation (Data S2A). However, last decades of the survey period (1984–1993 versus 2007–2016) nonlinear circular models showed that the rate at which flowering shows that the two lowest communities (1 and 2, n = 377 and 367 advanced accelerated over time for communities 1–4 (Data S2A). species, respectively) exhibited significant advances in flowering In all cases, nonlinear circular models (with year + year as phenology of 19.8 and 9.5 days, respectively (Figures 1A, 1B, predictors) fit better than linear models (with year as the only and 1F; Data S1B and S2C). In contrast, the two mid-elevation 434 Current Biology 30, 432–441, February 3, 2020 communities (3 and 4, n = 364 and 297 species, respectively) ex- amount of the variation in flowering time for low-, middle-, and hibited significant delays in flowering phenology of 7.4 and upper-elevation communities (1, 3, and 4); decreased precipita- 6.6 days, respectively (Figures 1C, 1D, and 1F; Data S1B and tion was associated with later flowering times (Table S1). For S2C). Within the highest community (5, n = 217 species), flower- communities 3 and 4, we detected significant interactions be- ing time significantly advanced by 1.9 days (Figures 1E and 1F; tween temperature and precipitation; flowering is expected to Data S1B and S2C). Delayed flowering times in communities 3 advance more when conditions are both warmer and wetter and 4 in the last decade in relation to the first, in combination (Table S1). The magnitudes of the predicted advances in flower- with the nonlinear trend toward earlier flowering over the full ing time associated with temperature decrease with increasing time series (Data S2A), indicate that flowering times in these elevation. For example, at the mean daily temperature for the communities were even later in the intervening years (1994– lowest community, an increase of 1 C is associated with a 2006) than in the most recent decade (Figure S1) and have 58.3-day advance in flowering phenology, whereas for the high- advanced at a relatively rapid rate in recent years. These shifts est community, an increase of 1 C is predicted to result in only a in flowering time translate into large differences in the extent to 2.25-day advance in flowering (Table S1). which mean flowering time differed among communities (Figures 2A, 2B, and 2C). For example, the interval between mean flower- Change in Flowering Synchrony of Subpopulations ing dates of the lowest- and the middle-elevation communities Our sample of subpopulations comprised 128 species, 77 of (1 and 3) increased by 27 days, almost an entire month, which occur in 4–5 elevation bands. At the subpopulation level, between the first and last decades (Data S2C). Flowering of com- we found significantly reduced synchrony of mean flowering munities at low to mid-elevations (2 and 3) became significantly over time for the 67 species in our sample found in both commu- less concentrated and more dispersed in time (community 2: nities 2 and 3. The difference in mean flowering times for subpop- 2 2 c = 9.64, p < 0.0019; community 3: c = 17.8, p < ulations in these two communities has become 3.26 days larger, 1 1 0.000025), as depicted by the shorter lengths of the arrows cor- growing from 2.52 to 5.66 days (Figure 4; Table S2). In contrast, responding to 2007–2016 in Figures 1B and 1C. synchrony significantly increased for subpopulations of the 30 species found in both communities 4 and 5, and mean flowering Change in Temperature and Precipitation became 3.90 days closer in time, shrinking from 4.88 to The rate of change in mean daily temperature differed among 1.05 days (Figure 4; Table S2). No changes in synchrony were elevations (significant interaction between year and elevation: detected for subpopulations found in communities 1 and 2 t = 2.53, p < 0.013), and temperature increased during the sur- (n = 56 species; Table S2) or 3 and 4 (n = 40 species; Table S2). vey period (1984–2016) at each elevation (971 m: R = 0.38, F = 18.7, p < 0.00015; 1,379 m: R = 0.49, F = 30.1, p < DISCUSSION 1,31 1,31 0.00001; 1,825 m: R = 0.58, F = 43.3, p < 0.00001; Figure 3A). 1,31 From the first decade (1984–1993) to the last (2007–2016), mean Although climate-change-driven shifts in flowering time have been daily temperature increased by 0.7 C at 957 m (from 19.6 C± widely documented, we lack understanding of the spatial variation 0.14 C [mean ± SE] to 20.3 C ± 0.11 C), by 1 C at 1,459 m in shifts within an ecosystem, leaving it unclear whether adjacent (from 17.0 C ± 0.18 C to 18.0 C ± 0.13 C), and by 1.3 Cat communities and subpopulations are shifting in unison or differ- 2,206 m (from 14.4 C ± 0.21 C to 15.7 C ± 0.14 C). Mean daily ently. By analyzing long-term phenological data that span an ele- temperature also increased over the longer term (1930–2016; vational gradient, we show that shifts in community-level flowering Figure S2A). In contrast, total annual precipitation decreased time are both spatially divergent and nonlinear over time. Across during the survey period, and the slopes did not vary among el- the entire 33-year time series, all but the highest community evations, although the intercepts did (conditional R = 0.51, shifted toward earlier flowering, but they did so at different rates. F = 21.1, p < 0.00001; Figure 3B). There was no significant Different rates of change occurred not only among different com- 1,96 linear trend in precipitation over the longer term (1930–2016; Fig- munities but also within communities, some showing accelerated ure S2B). However, from the first decade (1984–1993) to the last advances in flowering in more recent years. In the first versus the (2007–2016), total annual precipitation decreased by 26% at last decades of the time series, communities at lower elevations 957 m (from 44.4 ± 3.39 cm to 32.8 ± 1.93 cm), by 18% at shifted to flowering several weeks earlier, those at mid-elevations 1,459 m (from 60.0 ± 4.10 cm to 49.0 ± 3.75 cm), and by 26% shifted to flowering about a week later, and those at high elevation at 2,206 m (from 78.4 ± 5.06 cm to 58.3 ± 3.61 cm). shifted slightly earlier (Figure 1F). These results demonstrate that traditional calculations of longer-term trends, often made on the Changing Climate Drives Phenological Change basis of comparisons of two time points, can in fact mask more As described in the STAR Methods, climate variables were complex, nonlinear changes over time [41,42]. Within species, calculated for the 12-month window preceding each flowering synchrony in mean flowering phenology has decreased for sub- record. We did not detect temporal autocorrelation within the populations at mid-elevations while increasing for subpopulations precipitation or flowering-phenology time series (Ljung-Box at high elevations across the decades (Figure 4). Subpopulations tests: p > 0.05). Because temperature showed only weakly sig- can therefore differ significantly in their phenological responses nificant lag-1 autocorrelation (Ljung-Box test: p = 0.03), we did to changing climate conditions over small spatial scales [26,43], not perform any detrending prior to analyses. For all commu- potentially disrupting or augmenting gene flow and influencing nities, elevation-band-specific increases in temperature were local adaptation [27,28,44]. associated with earlier flowering times across the survey period Community-level flowering phenology becomes progressively (Table S1). Total 12-month precipitation explained a significant later with increasing elevation, but the difference in flowering Current Biology 30, 432–441, February 3, 2020 435 Figure 2. Change in Flowering Phenology of All Communities Comprising the Larger Metacommunity Circular plots show the dates of flowering, represented by colored lines, for all species in all communities for (A) 1984–1993 and (B) 2007–2016; the summary plot shows the change in the interval between mean flowering times for adjacent communities from 1984–1993 versus 2007–2016 (C). Within each circular plot, the lengths of the colored lines indicate the number of phenological observations on each date. Numbers give the day of year on which each month begins, and the arrows show the mean direction and mean resultant length; shorter arrows indicate more dispersed flowering times. (A) Orange, 945–1,079 m (community 1); yellow, 1,079–1,372 m (community 2); light green, 1,372–1,671 m (community 3); light blue, 1,671–1,939 m (community 4); pink, 1,939–2,212 m (community 5). (B) Red, 945–1,079 m (community 1); light orange, 1,079–1,372 m (community 2); dark green, 1,372–1,671 m (community 3); dark blue, 1,671–1,939 m (com- munity 4); purple, 1,939–2,212 m (community 5). See also Data S2A–S2C. time between the lower elevation bands grew ten days larger other trophic levels. Even without climate change, this semi- because communities at those elevations shifted at different arid environment constitutes a dynamic landscape of floral rates to earlier blooming (Figure 2; Data S2C). However, there resource availability that has a relatively low reliability of flower- was very little change in the interval between mean flowering ing [47]. Pollinators and other mobile species that forage across dates for the mid-elevation communities, and the difference in the elevational gradient could require behavioral adjustments to flowering time between the highest elevations actually became maintain temporal overlap with floral resources. The fact that 8 days smaller (Figure 2; Data S2C). Reduced differences in flowering shifted later between the first and last decades in the timing of leaf-out were also detected across elevations in some communities and earlier in others could mean that some the European Alps over six decades [45], and similarly spatially pollinators could extend their foraging seasons, which could complex shifts are occurring across latitudinal gradients [46]. have important implications for gene flow and reproductive Flowering times have also become more dispersed across the output of plants. calendar year for mid-elevation communities (Figures 1B and Whether species possess the ability to respond plastically to 1C). Together, these differential responses across space have ongoing climate change or must rely on adaptive evolution is important implications for patterns of resource availability for an open question [48–50]. Given overlap in species composition 436 Current Biology 30, 432–441, February 3, 2020 and precipitation tended to decrease with increasing elevation, as did the magnitude of observed shifts in flowering time. Thus, increasing temperatures and decreasing precipitation appear to be having less impact on communities at high elevation. By using all observations of flowering for each set of species and years, we avoided biases associated with the use of data on flowering onset alone [51]. Analyses that focus on flowering onset could be particularly problematic in ecosystems where flowering is sporadic or continuous. Although our analyses of subpopulation-level synchrony deal with changes in the mean flowering time, depending on the shape of the flowering curve, changes in the mean will reflect changes in peak flowering time and will likely influence interactions with mutualists and antagonists [52]. Ideally, however, analyses will encompass changes in the entire distributions of flowering, especially when the goal is to predict effects on gene flow. Even with data on the entire flowering record, the tendency for species to flower during both the latter and early parts of the calendar year mean that non-circular analyses will fail to accurately cap- ture the distribution of flowering times. Phenological events cross the calendar-year divide in many ecosystems, particularly those driven by precipitation in addition to temperature, such as the tropical, arid, and semi-arid ecosystems that cover most of the earth’s land area [53]. However, our current understanding of phenological responses to climate change is biased toward mid-latitude temperate ecosystems with discrete spring and summer flowering seasons [54]. Circular statistics as employed here provide powerful and underused alternative methods for analyzing phenological datasets [55] and yield novel insight into how phenological distributions are being affected by climate change. The environmental heterogeneity (sensu [56]) of the study area includes dissimilarities in land cover, vegetation, climate, hydrol- ogy, and topography both within and among communities. As a result, the mechanisms driving the spatially divergent and tempo- Figure 3. Change in Temperature and Precipitation rally nonlinear shifts in flowering phenology are likely to be many (A and B) (A) Significant positive relationships between mean daily temperature and interrelated but are difficult to determine given the paucity of and year and (B) significant negative relationships between total annual pre- cipitation and year, partitioned by elevation, for the survey period (1984–2016). research in semi-arid systems, particularly studies of large In (A), the slopes vary by elevation, whereas in (B), only the intercepts vary by numbers of species over elevational gradients. Observational elevation. Confidence intervals (95%) of regression slopes are shown. See also studies over elevational gradients are few [57], and experimental Figure S2 and Table S1. studies might not accurately predict plant phenological re- sponses [58]. Certainly, different species are likely to respond in among communities, with 28% of species spanning four or five diverse ways because of variation in sensitivity to changes in temperature and precipitation ([3]; ‘‘organismal mechanisms’’ elevation bands, our results suggest species can adjust flower- sensu [59]). Thus, the unique composition of species in each ing times in response to microclimates that vary across both space and time. In support of the idea that species are tracking community is likely partly responsible for divergent responses. elevation-specific microclimates, whether plastically, geneti- The topography of the study area and the microhabitats inhabited cally, or both, the modeled responses to yearly climate variables by the species in each community likely also influence the mech- align with the observed decadal differences in flowering time. For anisms involved [32]. Because of differences in exposure, evapo- example, the predicted advance of 58 days in mean flowering transpiration is likely greater for communities 1 and 2 than com- time per 1 C increase in temperature at the lowest elevation munities 3–5. Soils throughout the study area are uniformly band (Table S1) would generate an advance of 41 days with shallow lithosoils; organic matter, surface cover, nitrogen con- the 0.7 C increase in temperature that occurred between the first tent, and acidity tend to increase with elevation [60]. Additionally, and last decade (Figure 3A). This phenological advance is pre- coarse talus slopes holding pockets of deeper soils in commu- dicted to be countered by an expected 13-day delay in flowering nities 3 and 4 might retain more water than is possible in the lower generated by the 11.6-cm reduction in annual precipitation be- two communities, and the highly fractured bedrock could provide tween the decades (Figure 3B; Table S1), yielding a predicted water storage in community 5 if refreshed by precipitation [61]. net advance of 28 days, only 8 days more than the observed These differences in physical features and evapotranspiration advance. The predicted magnitude of the effects of temperature rates (‘‘environmental mechanisms’’ sensu [59]) could explain Current Biology 30, 432–441, February 3, 2020 437 Figure 4. Change in Flowering Synchrony of Subpopulations Change in the interval between mean flowering times for subpopulations within communities 1 versus 2, 2 versus 3, 3 versus 4, and 4 versus 5 during the first decade (1984–1993) and last decade (2007–2016) of the survey period. Positive values indicate longer intervals (decreased synchrony), and negative values indicate shorter intervals (increased synchrony). The difference in flowering time for subpopulations on elevation band 2 versus 3 was significantly larger in the last decade than in the first, whereas the difference in flowering time for subpopulations on elevation band 4 versus 5 was significantly smaller in the last decade than in the first. Highest posterior density intervals (95%) are shown. See also Table S2. changing climatic cues. Our study pro- vides a novel view of how the timing of flowering is changing in a semi-arid ecosystem, an ecosystem type that is ex- pected to expand with continued climate why flowering shifted at different rates across the gradient, such change [63]. Ephemeral and intermittent stream communities as communities 1 and 2 shifting earlier at rates 3–7 times faster in semi-arid ecosystems contain high biodiversity and provide than did communities 3 and 4. In addition, the responses of com- the same ecological services as do true riparian areas [64]; munities 3 and 4 were shaped by significant interactions between these systems, and the biodiversity they contain, are highly temperature and precipitation (Table S1), and flowering in those threatened by climate change [65]. The phenological changes communities advanced more when conditions were both warmer driven by increasing temperature and decreasing precipitation and wetter. Though correlational, this result suggests that these in the xeroriparian habitat in our study system could be indica- interacting drivers can give rise to divergent phenologies in neigh- tive of what we can expect elsewhere in the United States, and boring communities. particularly in the Southwest. This study shows that ecosystem Given the high variability of precipitation and recurring drought responses to climate change will be both variable and complex, in the Southwest [34], it is not surprising that we did not detect a particularly in highly heterogeneous systems characterized by long-term linear trend in total annual precipitation (Figure S2B). high interannual climatic variability, highly variable topography, However, since the late 1990s, the study area has been impacted and high biodiversity. Short-term studies of only a few by long-term warm hydraulic drought that is characterized by not species might not show the extent of change occurring, espe- only precipitation deficits but also rising temperatures [33,34]. cially when baseline data are lacking, and could in fact produce The considerable increase in precipitation variability in the South- erroneous conclusions. Our findings demonstrate that commu- west could exacerbate the effects of drought by reducing growth nities and subpopulations occupying different microclimates and increasing mortality [35]. In addition, the increasing intensity are exhibiting remarkably different responses to changing of monsoon storms [36] will likely result in greater run-off, which climatic conditions. The differences in both magnitude and di- means less moisture available to vegetation. In the study area, rection of responses highlight how climate change will result several dominant species have declined in abundance, particu- in community reshuffling in the temporal dimension. larly Carnegiea gigantea (saguaro), Juniperus deppeana var. dep- peana (alligator juniper), Parkinsonia microphylla (foothills palo STAR+METHODS verde), Pinus discolor (border pinyon), Pinus ponderosa subsp. brachyptera (ponderosa pine), and Quercus arizonica (Arizona Detailed methods are provided in the online version of this paper white oak), as have nearly all annuals, particularly when cool-sea- and include the following: son rains are poor [62]. The result of these declines is likely d KEY RESOURCES TABLE increased soil temperatures due to an increase in bare ground d LEAD CONTACT AND MATERIALS AVAILABILITY and decreased cover. These drought conditions could explain d EXPERIMENTAL MODEL AND SUBJECT DETAILS some of the temporally nonlinear responses we detected, leading d METHOD DETAILS communities 3 and 4 to exhibit delayed flowering in the last B Phenological Dataset decade of the time series in relation to the first. B Climate Data By virtue of long-term, taxonomically extensive, and highly d QUANTIFICATION AND STATISTICAL ANALYSIS temporally resolved data that span a spatial gradient, we B Circular Statistics were able to detect divergent metacommunity-level and sub- B Metacommunity Shifts in Flowering Phenology population-level shifts in flowering phenology driven by 438 Current Biology 30, 432–441, February 3, 2020 B Climate Models 11. Rafferty, N.E., Bertelsen, C.D., and Bronstein, J.L. (2016). Later flowering is associated with a compressed flowering season and reduced reproduc- B Subpopulation Changes in Flowering Synchrony tive output in an early season floral resource. Oikos 125, 821–828. d DATA AND CODE AVAILABILITY 12. Inouye, D.W. (2008). Effects of climate change on phenology, frost dam- age, and floral abundance of montane wildflowers. Ecology 89, 353–362. SUPPLEMENTAL INFORMATION 13. Iler, A.M., Compagnoni, A., Inouye, D.W., Williams, J.L., CaraDonna, P.J., Anderson, A., and Miller, T.E. (2019). Reproductive losses due to climate Supplemental Information can be found online at https://doi.org/10.1016/j. change-induced earlier flowering are not the primary threat to plant popu- cub.2019.11.071. lation viability in a perennial herb. J. Ecol. 107, 1931–1943. 14. Pearse, W.D., Davis, C.C., Inouye, D.W., Primack, R.B., and Davies, T.J. ACKNOWLEDGMENTS (2017). A statistical estimator for determining the limits of contemporary and historic phenology. Nat. Ecol. Evol. 1, 1876–1882. N.E.R. acknowledges UCR initial complement funding and thanks Paul Nabity for insightful discussions and David Rafferty for assistance with coding to 15. Kudo, G., and Cooper, E.J. (2019). When spring ephemerals fail to meet extract mean differences in flowering time within species. C.D.B. acknowl- pollinators: mechanism of phenological mismatch and its impact on plant edges the University of Arizona Herbarium staff who identified or verified nearly reproduction. Proc. Biol. Sci. 286, 20190573. all the taxa in the flora. We thank Joel Sachs and four anonymous reviewers for 16. Parsche, S., Frund, J., and Tscharntke, T. (2011). Experimental environ- comments that greatly improved the manuscript. mental change and mutualistic vs. antagonistic plant flower-visitor interac- tions. Perspect. Plant Ecol. Evol. Syst. 13, 27–35. AUTHOR CONTRIBUTIONS 17. Rafferty, N.E., and Ives, A.R. (2011). Effects of experimental shifts in flow- ering phenology on plant-pollinator interactions. Ecol. Lett. 14, 69–74. C.D.B. collected the data, N.E.R. analyzed the data with input from J.M.D., 18. Gezon, Z.J., Inouye, D.W., and Irwin, R.E. (2016). Phenological change in a N.E.R. wrote the first draft of the manuscript, and all authors devised the study spring ephemeral: implications for pollination and plant reproduction. questions and edited the manuscript. Glob. Change Biol. 22, 1779–1793. 19. Rafferty, N.E., and Ives, A.R. (2012). Pollinator effectiveness varies with DECLARATION OF INTERESTS experimental shifts in flowering time. Ecology 93, 803–814. 20. Forrest, J.R. (2015). Plant-pollinator interactions and phenological The authors declare no competing interests. change: what can we learn about climate impacts from experiments and observations? Oikos 124, 4–13. Received: August 25, 2019 Revised: October 20, 2019 21. Burkle, L.A., Marlin, J.C., and Knight, T.M. (2013). Plant-pollinator interac- Accepted: November 25, 2019 tions over 120 years: loss of species, co-occurrence, and function. Published: January 2, 2020 Science 339, 1611–1615. 22. Bartomeus, I., Ascher, J.S., Wagner, D., Danforth, B.N., Colla, S., REFERENCES Kornbluth, S., and Winfree, R. (2011). Climate-associated phenological advances in bee pollinators and bee-pollinated plants. Proc. Natl. Acad. 1. Bradley, N.L., Leopold, A.C., Ross, J., and Huffaker, W. (1999). Sci. USA 108, 20645–20649. Phenological changes reflect climate change in Wisconsin. Proc. Natl. 23. Thomson, J.D. (2019). Progressive deterioration of pollination service de- Acad. Sci. USA 96, 9701–9704. tected in a 17-year study vanishes in a 26-year study. New Phytol. 224, 2. Fitter, A.H., and Fitter, R.S. (2002). Rapid changes in flowering time in 1151–1159. British plants. Science 296, 1689–1691. 24. Yang, L.H., and Rudolf, V.H. (2010). Phenology, ontogeny and the effects 3. Crimmins, T.M., Crimmins, M.A., and Bertelsen, C.D. (2010). Complex re- of climate change on the timing of species interactions. Ecol. Lett. 13, sponses to climate drivers in onset of spring flowering across a semi-arid 1–10. elevation gradient. J. Ecol. 98, 1042–1051. 25. Aldridge, G., Inouye, D.W., Forrest, J.R., Barr, W.A., and Miller-Rushing, 4. Crimmins, T.M., Crimmins, M.A., and Bertelsen, C.D. (2011). Onset of A.J. (2011). Emergence of a mid-season period of low floral resources in summer flowering in a ‘Sky Island’ is driven by monsoon moisture. New a montane meadow ecosystem associated with climate change. J. Ecol. Phytol. 191, 468–479. 99, 905–913. 5. Calinger, K.M., Queenborough, S., and Curtis, P.S. (2013). Herbarium 26. Theobald, E.J., Breckheimer, I., and HilleRisLambers, J. (2017). Climate specimens reveal the footprint of climate change on flowering trends drives phenological reassembly of a mountain wildflower meadow com- across north-central North America. Ecol. Lett. 16, 1037–1044. munity. Ecology 98, 2799–2812. 6. Diez, J.M., Iba´ n˜ ez, I., Miller-Rushing, A.J., Mazer, S.J., Crimmins, T.M., 27. Ison, J.L., Wagenius, S., Reitz, D., and Ashley, M.V. (2014). Mating be- Crimmins, M.A., Bertelsen, C.D., and Inouye, D.W. (2012). Forecasting tween Echinacea angustifolia (Asteraceae) individuals increases with their phenology: from species variability to community patterns. Ecol. Lett. flowering synchrony and spatial proximity. Am. J. Bot. 101, 180–189. 15, 545–553. 28. Weis, A.E., Nardone, E., and Fox, G.A. (2014). The strength of assortative 7. Willis, C.G., Ruhfel, B., Primack, R.B., Miller-Rushing, A.J., and Davis, C.C. mating for flowering date and its basis in individual variation in flowering (2008). Phylogenetic patterns of species loss in Thoreau’s woods are schedule. J. Evol. Biol. 27, 2138–2151. driven by climate change. Proc. Natl. Acad. Sci. USA 105, 17029–17033. 29. Waser, N., and Price, M. (1991). Outcrossing distance effects in 8. Rafferty, N.E., and Nabity, P.D. (2017). A global test for phylogenetic signal Delphinium nelsonii: pollen loads, pollen tubes, and seed set. Ecology in shifts in flowering time under climate change. J. Ecol. 105, 627–633. 72, 171–179. 9. Miller-Rushing, A.J., Høye, T.T., Inouye, D.W., and Post, E. (2010). The ef- 30. Morellato, L., Alberton, B., Alvarado, S.T., Borges, B., Buisson, E., fects of phenological mismatches on demography. Philos. Trans. R. Soc. Camargo, M.G., Cancian, L.F., Carstensen, D.W., Escobar, D.F., Leite, Lond. B Biol. Sci. 365, 3177–3186. P.T., et al. (2016). Linking plant phenology to conservation biology. Biol. Conserv. 195, 60–72. 10. Cleland, E.E., Allen, J.M., Crimmins, T.M., Dunne, J.A., Pau, S., Travers, S.E., Zavaleta, E.S., and Wolkovich, E.M. (2012). Phenological tracking en- 31. Halbritter, A.H., Alexander, J.M., Edwards, P.J., and Billeter, R. (2013). ables positive species responses to climate change. Ecology 93, 1765– How comparable are species distributions along elevational and latitudinal 1771. climate gradients? Glob. Ecol. Biogeogr. 22, 1228–1237. Current Biology 30, 432–441, February 3, 2020 439 32. Ko¨ rner, C. (2007). The use of ‘altitude’ in ecological research. Trends Ecol. 49. Franks, S.J., Weber, J.J., and Aitken, S.N. (2014). Evolutionary and plastic Evol. 22, 569–574. responses to climate change in terrestrial plant populations. Evol. Appl. 7, 123–139. 33. Arizona State Climate Office (2019). Drought. Global Institute of 50. Merila, € J., and Hendry, A.P. (2014). Climate change, adaptation, and Sustainability, Arizona State University, Tempe, AZ. https://azclimate. phenotypic plasticity: the problem and the evidence. Evol. Appl. 7, 1–14. asu.edu/drought/. 51. Miller-Rushing, A.J., Inouye, D.W., and Primack, R.B. (2008). How well 34. Woodhouse, C.A., Meko, D.M., MacDonald, G.M., Stahle, D.W., and do first flowering dates measure plant responses to climate change? Cook, E.R. (2010). A 1,200-year perspective of 21st century drought in The effects of population size and sampling frequency. J. Ecol. 96, southwestern North America. Proc. Natl. Acad. Sci. USA 107, 21283– 1289–1296. 52. Elzinga, J.A., Atlan, A., Biere, A., Gigord, L., Weis, A.E., and Bernasconi, G. 35. Dannenberg, M.P., Wise, E.K., and Smith, W.K. (2019). Reduced tree (2007). Time after time: flowering phenology and biotic interactions. growth in the semiarid United States due to asymmetric responses to Trends Ecol. Evol. 22, 432–439. intensifying precipitation extremes. Sci. Adv. 5, w0667. 53. Kottek, M., Grieser, J., Beck, C., Rudolf, B., and Rubel, F. (2006). World 36. Demaria, E.M.C., Hazenberg, P., Scott, R.L., Meles, M.B., Nichols, M., and map of the Ko¨ ppen-Geiger climate classification updated. Meteorol. Z. Goodrich, D. (2019). Intensification of North American monsoon rainfall as (Berl.) 15, 259–263. observed from a long-term high-density gauge network. Geophys. Res. 54. Wolkovich, E.M., Cook, B.I., and Davies, T.J. (2014). Progress towards an Lett. 46, 6839–6847. interdisciplinary science of plant phenology: building predictions across 37. Crimmins, T.M., Crimmins, M.A., Bertelsen, D., and Balmat, J. (2008). space, time and species diversity. New Phytol. 201, 1156–1162. Relationships between alpha diversity of plant species in bloom and cli- 55. Morellato, L.P.C., Alberti, L.F., and Hudson, I.L. (2010). Applications of cir- matic variables across an elevation gradient. Int. J. Biometeorol. 52, cular statistics in plant phenology: a case studies approach. In 353–366. Phenological Research, L.L. Hudson, and M.R. Keatley, eds. (Dordrecht: 38. Crimmins, T.M., Crimmins, M., and Bertelsen, C.D. (2013). Temporal pat- Springer), pp. 339–359. terns in species flowering in Sky Islands of the Sonoran Desert ecoregion. 56. Stein, A., Gerstner, K., and Kreft, H. (2014). Environmental heterogeneity In Merging Science and Management in a Rapidly Changing World: as a universal driver of species richness across taxa, biomes and spatial Biodiversity and Management of the Madrean Archipelago III, G.J. scales. Ecol. Lett. 17, 866–880. Gottfried, P.F. Ffolliott, B.S. Gebow, and L.G. Eskew, eds. (U.S. 57. Sundqvist, M.K., Sanders, N.J., and Wardle, D.A. (2013). Community Department of Agriculture, Forest Service, Rocky Mountain Research and ecosystem responses to elevational gradients: processes, mecha- Station), pp. 33–39. nisms, and insights for global change. Annu. Rev. Ecol. Evol. Syst. 44, 39. Crimmins, T.M., Crimmins, M.A., and Bertelsen, C.D. (2013). Spring and 261–280. summer patterns in flowering onset, duration, and constancy across a wa- 58. Wolkovich, E.M., Cook, B.I., Allen, J.M., Crimmins, T.M., Betancourt, J.L., ter-limited gradient. Am. J. Bot. 100, 1137–1147. Travers, S.E., Pau, S., Regetz, J., Davies, T.J., Kraft, N.J., et al. (2012). 40. Crimmins, T.M., Bertelsen, D.C., and Crimmins, M.A. (2014). Within-sea- Warming experiments underpredict plant phenological responses to son flowering interruptions are common in the water-limited Sky Islands. climate change. Nature 485, 494–497. Int. J. Biometeorol. 58, 419–426. 59. Chmura, H.E., Kharouba, H.M., Ashander, J., Ehlman, S.M., Rivest, E.B., 41. Iler, A.M., Høye, T.T., Inouye, D.W., and Schmidt, N.M. (2013). Long-term and Yang, L.H. (2019). The mechanisms of phenology: the patterns and trends mask variation in the direction and magnitude of short-term pheno- processes of phenological shifts. Ecol. Monogr. 89, e01337–e22. logical shifts. Am. J. Bot. 100, 1398–1406. 60. Whittaker, R., Buol, S., Niering, W., and Havens, Y. (1968). A soil and vege- 42. Iler, A.M., Høye, T.T., Inouye, D.W., and Schmidt, N.M. (2013). Nonlinear tation pattern in the Santa Catalina Mountains, Arizona. Soil Sci. 105, flowering responses to climate: are species approaching their limits of 440–450. phenological change? Philos. Trans. R. Soc. Lond. B Biol. Sci. 368, 61. Schwinning, S. (2010). The ecohydrology of roots in rocks. Ecohydrology 3, 238–245. 43. de Keyzer, C.W., Rafferty, N.E., Inouye, D.W., and Thomson, J.D. (2017). 62. Bertelsen, C.D. (2018). Thirty-seven years on a mountain trail: vascular Confounding effects of spatial variation on shifts in phenology. Glob. flora and flowering phenology of the Finger Rock Canyon watershed, Change Biol. 23, 1783–1791. Santa Catalina Mountains, Arizona. Desert Plants 34, 3–268. 44. Hendry, A.P., and Day, T. (2005). Population structure attributable to 63. Rajaud, A., and de Noblet-Ducoudre,  N. (2017). Tropical semi-arid regions reproductive time: isolation by time and adaptation by time. Mol. Ecol. expanding over temperate latitudes under climate change. Clim. Change 14, 901–916. 144, 703–719. 45. Vitasse, Y., Signarbieux, C., and Fu, Y.H. (2018). Global warming leads to 64. Levick, L., Fonseca, J., Goodrich, D., Hernandez, M., Semmens, D., more uniform spring phenology across elevations. Proc. Natl. Acad. Sci. Stromberg, J., Leidy, R., Scianni, M., Guertin, D.P., Tluczek, M., et al. USA 115, 1004–1008. (2008). The ecological and hydrological significance of ephemeral and intermittent streams in the arid and semi-arid American Southwest. 46. Prevey, J., Vellend, M., Ru¨ ger, N., Hollister, R.D., Bjorkman, A.D., Myers- https://www.epa.gov/sites/production/files/2015-03/documents/ephemeral_ Smith, I.H., Elmendorf, S.C., Clark, K., Cooper, E.J., Elberling, B., et al. streams_report_final_508-kepner.pdf. (2017). Greater temperature sensitivity of plant phenology at colder sites: implications for convergence across northern latitudes. Glob. Change 65. Perry, L.G., Andersen, D.C., Reynolds, L.V., Nelson, S.M., and Shafroth, P.B. (2012). Vulnerability of riparian ecosystems to elevated CO2 and Biol. 23, 2660–2671. climate change in arid and semiarid western North America. Glob. 47. Wright, K.W., Vanderbilt, K.L., Inouye, D.W., Bertelsen, C.D., and Change Biol. 18, 821–842. Crimmins, T.M. (2015). Turnover and reliability of flower communities in 66. Brown, D.E. (1994). Biotic communities: Southwestern United States and extreme environments: insights from long-term phenology data sets. Northwestern Mexico (Salt Lake City: University of Utah Press). J. Arid Environ. 115, 27–34. 67. PRISM Climate Group Oregon State University (2017). created 25 Oct. 48. Anderson, J.T., Inouye, D.W., McKinney, A.M., Colautti, R.I., and Mitchell- http://prism.oregonstate.edu. Olds, T. (2012). Phenotypic plasticity and adaptive evolution contribute to advancing flowering phenology in response to climate change. Proc. Biol. 68. Daly, C., Halbleib, M., Smith, J.I., Gibson, W.P., Doggett, M.K., Taylor, Sci. 279, 3843–3852. G.H., Curtis, J., and Pasteris, P.P. (2008). Physiographically sensitive 440 Current Biology 30, 432–441, February 3, 2020 mapping of climatological temperature and precipitation across the 72. Cremers, J., Mulder, K.T., and Klugkist, I. (2018). Circular interpretation of conterminous United States. Int. J. Climatol. 28, 2031–2064. regression coefficients. Br. J. Math. Stat. Psychol. 71, 75–95. 69. Agostinelli, C., and Lund, U. (2017). R package ‘circular’: circular statistics 73. R Core Team (2019). R: A language and environment for statistical (version 0.4-93). URL. https://r-forge.r-project.org/projects/circular/. computing. R Foundation for Statistical Computing, Vienna, Austria. URL. https://www.R-project.org/. 70. Cremers, J. (2018). bpnreg: Bayesian projected normal regression models for circular data. R package version 1.0.0. URL. https://CRAN.R-project. 74. Pewsey, A., Neuhauser, M., and Ruxton, G.D. (2013). Circular Statistics in org/package=bpnreg. R (New York: Oxford University Press). 71. Cremers, J., and Klugkist, I. (2018). One direction? A tutorial for circular 75. Iler, A.M., Inouye, D.W., Schmidt, N.M., and Høye, T.T. (2017). Detrending data analysis using R with examples in cognitive psychology. Front. phenological time series improves climate-phenology analyses and re- Psychol. 9, 2040. veals evidence of plasticity. Ecology 98, 647–655. Current Biology 30, 432–441, February 3, 2020 441 STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Deposited Data Flowering phenology data This paper; Mendeley data Mendeley data (https://doi.org/10.17632/k6p34z78x9.1) Climate data This paper; Mendeley data Mendeley data (https://doi.org/10.17632/k6p34z78x9.1) LEAD CONTACT AND MATERIALS AVAILABILITY Further information and requests should be directed to and will be fulfilled by the Lead Contact, Nicole Rafferty (nicole.rafferty@ucr. edu). The datasets used in this study have been deposited to Mendeley data (https://doi.org/10.17632/k6p34z78x9.1). This study did not generate any new reagents. EXPERIMENTAL MODEL AND SUBJECT DETAILS The data come from systematic surveys by one coauthor (CDB) of all plant species and infraspecific taxa (hereafter ‘‘species’’) in flower along a trail in Finger Rock Canyon ascending to Mt. Kimball in the Santa Catalina Mountains of Arizona, USA (Figure S3). Although the canyon represents less than 1% of the area of the Santa Catalina range, 45% of the known plant taxa in the range have been found there [62]. In 8.05 km, the trail ascends from 945-2212 m, which was partitioned at the beginning of data collection by CDB into five elevation bands that include six vegetative associations in five biotic communities (based on [66]): 1) 945-1079 m (desert scrub, riparian scrub); 2) 1079-1372 m (desert scrub, scrub grassland); 3) 1372-1671 m (scrub grassland, oak woodland); 4) 1671-1939 m (oak-pine woodland); 5) 1939-2212 m (oak-pine woodland, pine forest; Figure S3). METHOD DETAILS Phenological Dataset Every species seen in anthesis (angiosperms) or releasing pollen (gymnosperms), together referred to as ‘‘flowering,’’ was recorded for each community along each 1.6-km-long trail segment on every survey. During the first nine years of data collection, a period characterized by above-normal precipitation, data were collected an average of 30 days per year, with at least two surveys per month during the growing seasons. Subsequently, data were collected an average of 50 days per year, nearly weekly. Because our analyses use all records of flowering and focus on mean flowering dates, this change in sampling frequency should not bias our estimates. Surveys were completed throughout the year with approximately 8% of the total number of surveys being completed each month of the year. The 33-year survey period (1984-2016) considered here comprises 169,030 observations collected during 1,639 surveys. Of the 590 species, 117 were observed in only one community, 140 were observed in two communities, 168 were observed in three communities, 100 were observed in four communities, and 65 were observed in all five communities. Additional details about the data collection protocol and transect can be found in Crimmins et al. [39] and Bertelsen [62]. In particular, Bertelsen [62] gives for each species the years flowering and the number of flowering observations per elevation band and month. Climate Data The primary source of climate data was the Parameter-elevation Regressions on Independent Slopes Model (PRISM) database [67], supplemented by on-site rain gauges. Gauges (Tru-Chek) were installed by one of us (CDB) in 2007 to obtain data specific to three locations: at 957 m (near the base of the transect), 1459 m (approximately midway up the transect), and 2206 m (near the peak of the transect). Each gauge was checked on average four times per month during 2007-2012 and 2014-2016, and mineral oil was used to prevent evaporation. Temperature patterns and precipitation data for years the gauge data were not available were extracted from 4 km PRISM cells that include the GPS coordinates of the gauges. PRISM data incorporates factors such as location, elevation, and topography in a climate-elevation regression for each grid cell [68]. Although two of the gauges are located within the same PRISM cell, GPS coordinates within the cell produce different values based on elevation. Monthly PRISM cell and rain gauge precipitation records are highly correlated for each of the three locations (r = 0.85-0.89). Thus, monthly temperature and precipitation data were extracted for the same approximate elevations as the gauge locations to create three elevation-specific climate predictors of flower- ing phenology at low-, mid-, and high-elevation. Based on detailed knowledge of the aspect and topography of each elevation band, long-term observation of weather patterns, and vegetation responses to short-term climatic events [see 62], we used low-elevation temperature and precipitation data to represent climate variables for communities 1 and 2, mid-elevation data to represent climate variables for community 3, and high-elevation data to represent climate variables for communities 4 and 5. Briefly, communities 1 and 2 have significant southern exposure, likely resulting in higher evapotranspiration, particularly since precipitation is less than at higher e1 Current Biology 30, 432–441.e1–e3, February 3, 2020 elevations. Community 3 is situated in the deepest and narrowest portion of the canyon, measured from the top of the ridges forming the canyon to the bottom of the drainage; temperatures are likely moderated by cold air drainage and the largely southwestern expo- sure. Community 4 has more continuous cover and a largely northwestern exposure, whereas community 5 has considerable cover but with extensive areas of exposed bedrock. The higher precipitation received in communities 3-5, the amount of cover, and expo- sures would likely lessen evapotranspiration. We regressed mean daily temperature and total annual precipitation against year to test whether these climate variables have changed over the survey period (1984-2016) and since 1930, when a sufficient number of nearby weather stations were available to provide reliable data for the study area [35]. In the initial regression models, we included the interaction between year and elevation to test whether slopes differed by elevation, in which case we fit separate regressions per elevation. If the interaction term was not significant, we fit a linear mixed-effects model with elevation as a random effect to allow intercepts to vary by elevation and compared the fit of models with and without this random effect (conditional versus marginal R ). QUANTIFICATION AND STATISTICAL ANALYSIS Circular Statistics Data from thirty years were used in the statistical analysis; 2004, 2005, and 2013 were excluded because surveys occurred irregularly during those years. We converted all survey dates to day of year (doy). For all circular statistics, we converted doy to radians and used the R packages ‘‘circular’’ version 0.4-93 [69] and, to construct circular mixed-effects models, ‘‘bpnreg’’ version 1.0.0 [70]. Additional details regarding circular mixed-effects models, their formulation, and interpretation can be found in Cremers and Klugkist [71] and Cremers et al. [72]. Briefly, for the circular mixed-effects models, statistical significance of continuous predictors was gauged by whether the 95% highest posterior density (HPD) lower and upper bounds for circular model coefficients included zero (not signif- icant); significance of categorical variables was judged by whether the 95% HPD intervals for both component I and component II linear coefficients included zero (not significant; [72]). For circular differences between variables, significance was determined by the proportion of iterations that were negative. Models were compared using the deviance information criterion (DIC and DIC ) alt and the Watanabe-Akaike information criterion (WAIC and WAIC ), which both reward better-fitting models while penalizing model 1 2 complexity [70]. All analyses were conducted in R version 3.5.2 [73]. Metacommunity Shifts in Flowering Phenology To examine how flowering phenology at the metacommunity level has changed over time while holding elevation constant, we compared all dates on which flowers of any species were observed within a given community from 1984-2016. This analysis allowed us to determine whether communities within each elevation band exhibited trends of advanced or delayed flowering across all years. We fitted additional models that included a quadratic term for year to test for nonlinear changes in flowering phenology over the full time series. We were also interested in examining trends in the early versus later years of the survey period with the expectation that any phenological changes in response to climate change would be most apparent when comparing the two end segments of the time series. In particular, phenological effects of the drought that began in the late 1990s are likely to be apparent in the latest years of the time series. Therefore, we also examined how flowering time changed per community between the first and last decades of the data- set (i.e., 1984-1993 versus 2007-2016), which additionally enabled us to visualize any shifts with circular plots. Because these ana- lyses use all available flowering records for each community, the various flowering distributions (including any bimodal distributions) are aggregated when mean flowering times are calculated. For each community, we constructed circular mixed-effects models with doy of flowering (in radians) as the response and year (continuous), year + year , or decade (categorical) as the predictor(s), with spe- cies identity included as a random effect to account for repeated observations of the same species over time. Year was centered at zero to aid interpretation of model coefficients. We also verified that communities exhibited elevation-band-specific responses in phenology by fitting a circular mixed-effects model with year, elevation band (continuous), and the interaction between year and elevation band as predictors (and species identity as a random effect). The bpnme () function within the bpnreg package uses a Bayesian approach and Markov chain Monte Carlo (MCMC) samplers to estimate model parameters [71]. For each model, we ran 10,000 iterations with a burn-in period of 1,000 iterations and no lag because there was minimal auto-correlation detected in the MCMC chains. Because the model with the interaction between year and elevation band as a predictor was very computationally intensive, we ran only 1,000 iterations with a burn-in period of 100 iterations for that model. We inspected traceplots to verify that models had converged, which is the method recommended by the package authors [70]. For models with continuous predictors, we report the ‘slope at the mean’ (SAM) circular model coefficients because they are the least biased [72]. To test whether the con- centration of flowering times (i.e., the spread of flowering times throughout the calendar year) changed for a given community be- tween decades, we used Wallraff’s test for a common concentration [74]. Climate Models To test whether climate variables were related to observed shifts in flowering time within each community, we used circular mixed- effects models with doy of flowering as the response and mean daily temperature and total precipitation during the 12-month window preceding and including the month in which flowering was observed as predictors, with species identity as a random effect. In our initial models, we also included the interaction between temperature and precipitation. These climate variables are specific to each flowering record and capture temperature and precipitation conditions for a one-year period before each observation, regardless of the calendar date of flowering. To test for temporal autocorrelation within each time series (precipitation, temperature, and flowering Current Biology 30, 432–441.e1–e3, February 3, 2020 e2 phenology), we used Ljung-Box tests with a lag of one year [75]. Temperature and precipitation variables were centered at zero to aid interpretation of model coefficients. As before, we ran 10,000 iterations per model with a burn-in period of 1,000 iterations and no lag because there was minimal auto-correlation detected in the MCMC chains. We inspected traceplots to verify that MCMC chains had converged [70]. For these models, we report the SAM circular model coefficients [72]. Subpopulation Changes in Flowering Synchrony To analyze within-species shifts in flowering phenology across time and space, we limited the dataset to only those species that: (i) occurred in at least two adjacent communities along the transect (i.e., communities 1 and 2; communities 2 and 3; communities 3 and 4; communities 4 and 5), (ii) had been observed flowering during at least four years in the first decade (1984-1993) and four years in the last decade (2007-2016), and (iii) had been observed at least four times per community in each year flowering was documented. The resulting dataset comprised 128 unique species with subpopulations in two or more adjacent communities. We then calculated the mean doy of flowering (in radians) per species per community per decade and took the difference between these means for each community per decade (e.g., we subtracted the mean doy of flowering for a given species in community 2 in the first decade from the mean doy of flowering for the same species in community 1 in the first decade; we repeated this process for the last decade). These values indicate how much the mean flowering times of subpopulations differed in the first versus the last decade of the dataset and provide a way to test for changes in subpopulation-level synchrony. These differences were first converted back to a circular variable so that values corresponded to the number of days in radians on the unit circle, measuring counterclockwise from 0 radians if the differences were positive and clockwise from 2p radians if the differences were negative. We then used these differences as our response variable in circular mixed-effects models with decade as the predictor and species identity as a random effect. We used the same model specifications as in previous models. DATA AND CODE AVAILABILITY The phenology and climate datasets used in this study are available at Mendeley data (https://doi.org/10.17632/k6p34z78x9.1). e3 Current Biology 30, 432–441.e1–e3, February 3, 2020 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Current Biology Unpaywall

Changing Climate Drives Divergent and Nonlinear Shifts in Flowering Phenology across Elevations

Current BiologyFeb 1, 2020

Loading next page...
 
/lp/unpaywall/changing-climate-drives-divergent-and-nonlinear-shifts-in-flowering-6P2shvt5LO

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Unpaywall
ISSN
0960-9822
DOI
10.1016/j.cub.2019.11.071
Publisher site
See Article on Publisher Site

Abstract

Article Changing Climate Drives Divergent and Nonlinear Shifts in Flowering Phenology across Elevations Highlights Authors d Flowering time is diverging among communities across an Nicole E. Rafferty, Jeffrey M. Diez, elevational gradient C. David Bertelsen d Divergence reflects nonlinear shifts in flowering phenology Correspondence over three decades nicole.rafferty@ucr.edu d Climatic variables have also changed differently across In Brief space and over time Using 33 years of data, Rafferty et al. test d Changing climate is driving phenological reshuffling across whether flowering times have diverged local spatial gradients for neighboring communities comprising 590 species across an elevational gradient. Divergent and nonlinear shifts in flowering time are related to elevation- specific changes in temperature and precipitation and will likely alter species interactions and gene flow. Rafferty et al., 2020, Current Biology 30, 432–441 February 3, 2020 ª 2019 The Author(s). Published by Elsevier Ltd. https://doi.org/10.1016/j.cub.2019.11.071 Current Biology Article Changing Climate Drives Divergent and Nonlinear Shifts in Flowering Phenology across Elevations 1,2,6, 3 4,5 Nicole E. Rafferty, * Jeffrey M. Diez, and C. David Bertelsen Department of Evolution, Ecology, and Organismal Biology, University of California, Riverside, 900 University Avenue, Riverside, CA 92521, USA Rocky Mountain Biological Laboratory, PO Box 519, Crested Butte, CO 81224, USA Department of Botany and Plant Sciences, University of California, Riverside, 900 University Avenue, Riverside, CA 92521, USA School of Natural Resources and the Environment, University of Arizona, 1955 E. Sixth Street, Tucson, AZ 85721, USA Herbarium, University of Arizona, PO Box 210036, Tucson, AZ 85721, USA Lead Contact *Correspondence: nicole.rafferty@ucr.edu https://doi.org/10.1016/j.cub.2019.11.071 SUMMARY climate is driving phenological reshuffling across local spatial gradients. Climate change is known to affect regional weather patterns and phenology; however, we lack under- INTRODUCTION standing of how climate drives phenological change across local spatial gradients. This spatial variation is Flowering phenology is a key biological indicator of current climate change [1–5]. However, the magnitude and direction of critical for determining whether subpopulations changes in flowering times appear to vary significantly both and metacommunities are changing in unison or among species and among communities, in part because of vari- diverging in phenology. Divergent responses could ation in how much abiotic factors such as temperature and pre- reduce synchrony both within species (disrupting cipitation have changed and in how sensitive species are to gene flow among subpopulations) and among spe- those changes [6]. Variation in responses has also been linked cies (disrupting interspecific interactions in commu- to biotic factors. For example, in some communities, flowering nities). We also lack understanding of phenological phenology responses are related to plant traits, such as flower- change in environments where life history events ing season [5], whether species have an annual or perennial life are frequently aseasonal, such as the tropical, arid, cycle [2,5], and whether species are wind or animal pollinated and semi-arid ecosystems that cover vast areas. [5]. In general, species that flower early in the season or are an- Using a 33-year-long dataset spanning a 1,267-m nuals show the greatest advances in phenology, whereas evi- dence regarding pollination mode is mixed [2,5]. In some cases, semi-arid elevational gradient in the southwestern closely related species exhibit similar responses; thus, phyloge- United States, we test whether flowering phenology netic relationships can also be important predictors of shifts in diverged among subpopulations within species flowering time [7,8]. and among five communities comprising 590 spe- Altered flowering phenology in response to climate can affect cies. Applying circular statistics to test for changes population dynamics and demography via several avenues [9]. in year-round flowering, we show flowering has For example, plant species that track climate by advancing become earlier for all communities except at the phenologically have higher metrics of performance such as highest elevations. However, flowering times shifted individual growth [10]. Delayed flowering has been linked to at different rates across elevations likely because compression of the flowering period and lower fruit and seed of elevation-specific changes in temperature and set [11]. The timing of flowering in relation to snowmelt and precipitation, indicating diverging phenologies of damaging frost events can have large effects on floral abun- dance in natural populations of wildflowers [12], although conse- neighboring communities. Subpopulations of indi- quent floral abortion might not translate into reduced population vidual species also diverged at mid-elevation but viability [13]. Nonetheless, phenological responsiveness has converged in phenology at high elevation. These been shown to be an important target of natural selection [8]. changes in flowering phenology among communities Increased variability in flowering phenology among species and subpopulations are undetectable when data are over time has been detected in some temperate datasets [14] pooled across the gradient. Furthermore, we show and is also likely to have demographic consequences. For that nonlinear changes in flowering times over the example, greater variation in species’ flowering times could alter 33-year record are obscured by traditional calcula- temporal overlap among different species, potentially affecting tions of long-term trends. These findings reveal pollen transfer and seed set. greater spatiotemporal complexity in phenological While communities undergo temporal reshuffling, shifts responses than previously recognized and indicate in flowering phenology can alter interspecific interactions, as 432 Current Biology 30, 432–441, February 3, 2020 ª 2019 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). well. For example, earlier flowering can cause asynchrony be- variation in temperature [31], although key abiotic factors vary tween plants and pollinators, particularly for species that flower in different ways along elevational and latitudinal gradients in early spring and are visited by the earliest emerging pollina- [32]. Our primary goal is to determine whether flowering phenol- tors, such as queen bumble bees [15]. Experimental manipula- ogies of communities and subpopulations are shifting in the tions of phenology have demonstrated that flowering time can same ways across elevations or are responding differently over affect the frequency of both plant-herbivore [16] and plant-polli- space. Differential shifts at the community or subpopulation nator interactions [17,18] and how effectively plants are polli- levels would signal changes in temporal overlap or synchrony nated ([19], but also see [20]). Over longer timescales, phenolog- across elevations. A second goal is to determine whether chang- ical mismatches between plants and pollinators can account for ing temperature and precipitation patterns are responsible for the loss of some interactions and have likely contributed to local any community-specific changes in phenology. Since the late species extinctions [21], although in other cases, synchrony be- 1990s, the study area has been impacted by a long-term warm tween plants and pollinators is maintained under climate change drought that is characterized by not only precipitation deficits [22]. Indeed, analyses of pollination success over very long time- but also rising temperatures, increasing precipitation intensity, scales could be necessary to determine trends [23]. Competitive and increasing precipitation variability [33–36]. This drought and facilitative interactions between plant species are also likely and longer-term climatic changes could be affecting commu- to be reshaped by phenological shifts [24]. nities and subpopulations along this gradient differently, in part Overall, research to date has provided some understanding of because of species- and population-specific responses and what traits predict species-level flowering phenology shifts and spatial environmental variation that could modulate the changing how those shifts influence performance, demography, and spe- abiotic conditions. Thus, we seek to determine whether eleva- cies interactions. However, how long-term shifts in flowering tional differences in climate variables and their changes over phenology vary among species-rich communities and subpopu- time lead to temporally nonlinear shifts in community-level lations across semi-arid environmental gradients has not been phenology that could result in spatially divergent phenological previously investigated. Prior work on change in flowering changes. phenology across spatial gradients has instead addressed The dataset we analyze is unique; in addition to capturing a topics such as how shifts in flowering time vary among species steep elevational gradient comprising five communities (defined [3], across a mosaic of moisture habitats at a single elevation by elevation bands) with six intergrading vegetative associa- [25], or over only a few years [26]. Ultimately, divergence in flow- tions, it includes 590 vascular plant taxa ranging from annuals ering phenology among subpopulations (i.e., spatially structured to trees and was collected in a semi-arid ecosystem where pre- subsets of larger populations) is expected to result in decreased cipitation and temperature are key triggers of flowering for many pollen flow and greater reproductive isolation because conspe- species at low and high elevation, respectively [3,4,37]. There is cifics overlap less in flowering time across space [27,28]. little overlap in the flowering assemblages of the five commu- Equally, convergence in phenology via increased overlap in flow- nities except in winter [37]. Fifty percent or more of species in ering time across space could result in increased pollen flow and each community are highly opportunistic, wherein flowering is outbreeding depression if subpopulations are locally adapted triggered by antecedent climatic conditions [38–40]. Because [29]. At the metacommunity level, divergent phenological re- these species have one to several flowering periods of varying sponses among adjacent communities will likely affect trophic durations, and flowering late in one calendar year could be the and non-trophic interactions, leading to altered community early part of a flowering season continuing into the next calendar structure and ecosystem processes [30]. For example, if flower- year, we apply circular statistical methods to detect spatiotem- ing phenology in a montane plant community shifts unevenly poral patterns. across elevations, this will alter the timing of resource availability We address the following questions: (1) at the metacommunity for species in other trophic levels, such as pollinators, and will level, has flowering phenology diverged or converged across el- likely have downstream effects on seed and fruit production, evations; (2) can variation in temperature and precipitation affecting frugivores, plant recruitment, and competitive interac- across elevations explain shifts in flowering phenology; and, (3) tions. Despite the many possible ways in which climate- within species, has flowering become more, or less, synchronous change-driven shifts in flowering phenology could affect sub- for subpopulations across space? Our approach reveals spatially populations and metacommunities, investigation of these divergent, temporally nonlinear changes in flowering phenology topics is limited by a lack of long-term phenological data across among communities and subpopulations that are obscured by the spatial scales important for maintaining gene flow and the pooling of phenological data across space and the simple species interactions. Without such data, it is impossible to deter- linear calculation of long-term trends. These results highlight mine whether species are responding at finer spatial scales that the utility of circular statistics for detecting patterns in phenology could alter the temporal overlap among subpopulations and and point to a greater complexity of responses than has previ- communities. ously been recognized, suggesting that climate change will Here, we analyze a 33-year record of flowering phenology that reshuffle communities in multiple dimensions. spans an elevational gradient of 1,267 m in the southwestern United States that encompasses desert scrub at low elevation, RESULTS several semi-arid associations (riparian scrub, scrub grassland, oak woodland, and oak-pine woodland), and pine forest at Change in Flowering Phenology from 1984–2016 high elevation. This elevational gradient captures a range of The relationship between flowering time and year varied environmental conditions that could be indicative of latitudinal by elevation band, indicating elevation-specific phenological Current Biology 30, 432–441, February 3, 2020 433 Figure 1. Change in Flowering Phenology of Each Community (A–E) Circular plots show the dates of flowering, represented by colored lines, for two time periods for all species within each community ascending the transect. Within each circular plot, the lengths of the colored lines indicate the number of phenological observations on each date. The lighter color represents data from 1984–1993, and the darker color represents data from 2007–2016; numbers give the day of year on which each month begins. The arrows show the mean di- rection and the mean resultant length (a metric of concentration) for a given time period; shorter arrows indicate more dispersed flowering times. (A) 945–1,079 m (community 1); (B) 1,079–1,372 m (community 2); (C) 1,372–1,671 m (community 3); (D) 1,671–1,939 m (community 4); (E) 1,939–2,212 m (community 5). (F) The summary plot shows the change in flowering time for each community. See also Figure S1; Data S1A and S1B. responses of communities (Data S1A). We therefore report re- predictor; Data S2B). No significant change in flowering time sults per community, numbering communities from 1–5 and was detected for the highest-elevation community (5; Data S2A). ordering from the lowest (1) to highest (5) elevation bands. When examined with linear circular models, from 1984–2016, Change in Flowering Phenology from 1984–1993 versus communities 1–4 shifted significantly earlier at a rate of 2.5, 2007–2016 1.6, 0.44, and 0.36 days per year, respectively, and shifts An examination of changes in flowering time during the first and became smaller with increasing elevation (Data S2A). However, last decades of the survey period (1984–1993 versus 2007–2016) nonlinear circular models showed that the rate at which flowering shows that the two lowest communities (1 and 2, n = 377 and 367 advanced accelerated over time for communities 1–4 (Data S2A). species, respectively) exhibited significant advances in flowering In all cases, nonlinear circular models (with year + year as phenology of 19.8 and 9.5 days, respectively (Figures 1A, 1B, predictors) fit better than linear models (with year as the only and 1F; Data S1B and S2C). In contrast, the two mid-elevation 434 Current Biology 30, 432–441, February 3, 2020 communities (3 and 4, n = 364 and 297 species, respectively) ex- amount of the variation in flowering time for low-, middle-, and hibited significant delays in flowering phenology of 7.4 and upper-elevation communities (1, 3, and 4); decreased precipita- 6.6 days, respectively (Figures 1C, 1D, and 1F; Data S1B and tion was associated with later flowering times (Table S1). For S2C). Within the highest community (5, n = 217 species), flower- communities 3 and 4, we detected significant interactions be- ing time significantly advanced by 1.9 days (Figures 1E and 1F; tween temperature and precipitation; flowering is expected to Data S1B and S2C). Delayed flowering times in communities 3 advance more when conditions are both warmer and wetter and 4 in the last decade in relation to the first, in combination (Table S1). The magnitudes of the predicted advances in flower- with the nonlinear trend toward earlier flowering over the full ing time associated with temperature decrease with increasing time series (Data S2A), indicate that flowering times in these elevation. For example, at the mean daily temperature for the communities were even later in the intervening years (1994– lowest community, an increase of 1 C is associated with a 2006) than in the most recent decade (Figure S1) and have 58.3-day advance in flowering phenology, whereas for the high- advanced at a relatively rapid rate in recent years. These shifts est community, an increase of 1 C is predicted to result in only a in flowering time translate into large differences in the extent to 2.25-day advance in flowering (Table S1). which mean flowering time differed among communities (Figures 2A, 2B, and 2C). For example, the interval between mean flower- Change in Flowering Synchrony of Subpopulations ing dates of the lowest- and the middle-elevation communities Our sample of subpopulations comprised 128 species, 77 of (1 and 3) increased by 27 days, almost an entire month, which occur in 4–5 elevation bands. At the subpopulation level, between the first and last decades (Data S2C). Flowering of com- we found significantly reduced synchrony of mean flowering munities at low to mid-elevations (2 and 3) became significantly over time for the 67 species in our sample found in both commu- less concentrated and more dispersed in time (community 2: nities 2 and 3. The difference in mean flowering times for subpop- 2 2 c = 9.64, p < 0.0019; community 3: c = 17.8, p < ulations in these two communities has become 3.26 days larger, 1 1 0.000025), as depicted by the shorter lengths of the arrows cor- growing from 2.52 to 5.66 days (Figure 4; Table S2). In contrast, responding to 2007–2016 in Figures 1B and 1C. synchrony significantly increased for subpopulations of the 30 species found in both communities 4 and 5, and mean flowering Change in Temperature and Precipitation became 3.90 days closer in time, shrinking from 4.88 to The rate of change in mean daily temperature differed among 1.05 days (Figure 4; Table S2). No changes in synchrony were elevations (significant interaction between year and elevation: detected for subpopulations found in communities 1 and 2 t = 2.53, p < 0.013), and temperature increased during the sur- (n = 56 species; Table S2) or 3 and 4 (n = 40 species; Table S2). vey period (1984–2016) at each elevation (971 m: R = 0.38, F = 18.7, p < 0.00015; 1,379 m: R = 0.49, F = 30.1, p < DISCUSSION 1,31 1,31 0.00001; 1,825 m: R = 0.58, F = 43.3, p < 0.00001; Figure 3A). 1,31 From the first decade (1984–1993) to the last (2007–2016), mean Although climate-change-driven shifts in flowering time have been daily temperature increased by 0.7 C at 957 m (from 19.6 C± widely documented, we lack understanding of the spatial variation 0.14 C [mean ± SE] to 20.3 C ± 0.11 C), by 1 C at 1,459 m in shifts within an ecosystem, leaving it unclear whether adjacent (from 17.0 C ± 0.18 C to 18.0 C ± 0.13 C), and by 1.3 Cat communities and subpopulations are shifting in unison or differ- 2,206 m (from 14.4 C ± 0.21 C to 15.7 C ± 0.14 C). Mean daily ently. By analyzing long-term phenological data that span an ele- temperature also increased over the longer term (1930–2016; vational gradient, we show that shifts in community-level flowering Figure S2A). In contrast, total annual precipitation decreased time are both spatially divergent and nonlinear over time. Across during the survey period, and the slopes did not vary among el- the entire 33-year time series, all but the highest community evations, although the intercepts did (conditional R = 0.51, shifted toward earlier flowering, but they did so at different rates. F = 21.1, p < 0.00001; Figure 3B). There was no significant Different rates of change occurred not only among different com- 1,96 linear trend in precipitation over the longer term (1930–2016; Fig- munities but also within communities, some showing accelerated ure S2B). However, from the first decade (1984–1993) to the last advances in flowering in more recent years. In the first versus the (2007–2016), total annual precipitation decreased by 26% at last decades of the time series, communities at lower elevations 957 m (from 44.4 ± 3.39 cm to 32.8 ± 1.93 cm), by 18% at shifted to flowering several weeks earlier, those at mid-elevations 1,459 m (from 60.0 ± 4.10 cm to 49.0 ± 3.75 cm), and by 26% shifted to flowering about a week later, and those at high elevation at 2,206 m (from 78.4 ± 5.06 cm to 58.3 ± 3.61 cm). shifted slightly earlier (Figure 1F). These results demonstrate that traditional calculations of longer-term trends, often made on the Changing Climate Drives Phenological Change basis of comparisons of two time points, can in fact mask more As described in the STAR Methods, climate variables were complex, nonlinear changes over time [41,42]. Within species, calculated for the 12-month window preceding each flowering synchrony in mean flowering phenology has decreased for sub- record. We did not detect temporal autocorrelation within the populations at mid-elevations while increasing for subpopulations precipitation or flowering-phenology time series (Ljung-Box at high elevations across the decades (Figure 4). Subpopulations tests: p > 0.05). Because temperature showed only weakly sig- can therefore differ significantly in their phenological responses nificant lag-1 autocorrelation (Ljung-Box test: p = 0.03), we did to changing climate conditions over small spatial scales [26,43], not perform any detrending prior to analyses. For all commu- potentially disrupting or augmenting gene flow and influencing nities, elevation-band-specific increases in temperature were local adaptation [27,28,44]. associated with earlier flowering times across the survey period Community-level flowering phenology becomes progressively (Table S1). Total 12-month precipitation explained a significant later with increasing elevation, but the difference in flowering Current Biology 30, 432–441, February 3, 2020 435 Figure 2. Change in Flowering Phenology of All Communities Comprising the Larger Metacommunity Circular plots show the dates of flowering, represented by colored lines, for all species in all communities for (A) 1984–1993 and (B) 2007–2016; the summary plot shows the change in the interval between mean flowering times for adjacent communities from 1984–1993 versus 2007–2016 (C). Within each circular plot, the lengths of the colored lines indicate the number of phenological observations on each date. Numbers give the day of year on which each month begins, and the arrows show the mean direction and mean resultant length; shorter arrows indicate more dispersed flowering times. (A) Orange, 945–1,079 m (community 1); yellow, 1,079–1,372 m (community 2); light green, 1,372–1,671 m (community 3); light blue, 1,671–1,939 m (community 4); pink, 1,939–2,212 m (community 5). (B) Red, 945–1,079 m (community 1); light orange, 1,079–1,372 m (community 2); dark green, 1,372–1,671 m (community 3); dark blue, 1,671–1,939 m (com- munity 4); purple, 1,939–2,212 m (community 5). See also Data S2A–S2C. time between the lower elevation bands grew ten days larger other trophic levels. Even without climate change, this semi- because communities at those elevations shifted at different arid environment constitutes a dynamic landscape of floral rates to earlier blooming (Figure 2; Data S2C). However, there resource availability that has a relatively low reliability of flower- was very little change in the interval between mean flowering ing [47]. Pollinators and other mobile species that forage across dates for the mid-elevation communities, and the difference in the elevational gradient could require behavioral adjustments to flowering time between the highest elevations actually became maintain temporal overlap with floral resources. The fact that 8 days smaller (Figure 2; Data S2C). Reduced differences in flowering shifted later between the first and last decades in the timing of leaf-out were also detected across elevations in some communities and earlier in others could mean that some the European Alps over six decades [45], and similarly spatially pollinators could extend their foraging seasons, which could complex shifts are occurring across latitudinal gradients [46]. have important implications for gene flow and reproductive Flowering times have also become more dispersed across the output of plants. calendar year for mid-elevation communities (Figures 1B and Whether species possess the ability to respond plastically to 1C). Together, these differential responses across space have ongoing climate change or must rely on adaptive evolution is important implications for patterns of resource availability for an open question [48–50]. Given overlap in species composition 436 Current Biology 30, 432–441, February 3, 2020 and precipitation tended to decrease with increasing elevation, as did the magnitude of observed shifts in flowering time. Thus, increasing temperatures and decreasing precipitation appear to be having less impact on communities at high elevation. By using all observations of flowering for each set of species and years, we avoided biases associated with the use of data on flowering onset alone [51]. Analyses that focus on flowering onset could be particularly problematic in ecosystems where flowering is sporadic or continuous. Although our analyses of subpopulation-level synchrony deal with changes in the mean flowering time, depending on the shape of the flowering curve, changes in the mean will reflect changes in peak flowering time and will likely influence interactions with mutualists and antagonists [52]. Ideally, however, analyses will encompass changes in the entire distributions of flowering, especially when the goal is to predict effects on gene flow. Even with data on the entire flowering record, the tendency for species to flower during both the latter and early parts of the calendar year mean that non-circular analyses will fail to accurately cap- ture the distribution of flowering times. Phenological events cross the calendar-year divide in many ecosystems, particularly those driven by precipitation in addition to temperature, such as the tropical, arid, and semi-arid ecosystems that cover most of the earth’s land area [53]. However, our current understanding of phenological responses to climate change is biased toward mid-latitude temperate ecosystems with discrete spring and summer flowering seasons [54]. Circular statistics as employed here provide powerful and underused alternative methods for analyzing phenological datasets [55] and yield novel insight into how phenological distributions are being affected by climate change. The environmental heterogeneity (sensu [56]) of the study area includes dissimilarities in land cover, vegetation, climate, hydrol- ogy, and topography both within and among communities. As a result, the mechanisms driving the spatially divergent and tempo- Figure 3. Change in Temperature and Precipitation rally nonlinear shifts in flowering phenology are likely to be many (A and B) (A) Significant positive relationships between mean daily temperature and interrelated but are difficult to determine given the paucity of and year and (B) significant negative relationships between total annual pre- cipitation and year, partitioned by elevation, for the survey period (1984–2016). research in semi-arid systems, particularly studies of large In (A), the slopes vary by elevation, whereas in (B), only the intercepts vary by numbers of species over elevational gradients. Observational elevation. Confidence intervals (95%) of regression slopes are shown. See also studies over elevational gradients are few [57], and experimental Figure S2 and Table S1. studies might not accurately predict plant phenological re- sponses [58]. Certainly, different species are likely to respond in among communities, with 28% of species spanning four or five diverse ways because of variation in sensitivity to changes in temperature and precipitation ([3]; ‘‘organismal mechanisms’’ elevation bands, our results suggest species can adjust flower- sensu [59]). Thus, the unique composition of species in each ing times in response to microclimates that vary across both space and time. In support of the idea that species are tracking community is likely partly responsible for divergent responses. elevation-specific microclimates, whether plastically, geneti- The topography of the study area and the microhabitats inhabited cally, or both, the modeled responses to yearly climate variables by the species in each community likely also influence the mech- align with the observed decadal differences in flowering time. For anisms involved [32]. Because of differences in exposure, evapo- example, the predicted advance of 58 days in mean flowering transpiration is likely greater for communities 1 and 2 than com- time per 1 C increase in temperature at the lowest elevation munities 3–5. Soils throughout the study area are uniformly band (Table S1) would generate an advance of 41 days with shallow lithosoils; organic matter, surface cover, nitrogen con- the 0.7 C increase in temperature that occurred between the first tent, and acidity tend to increase with elevation [60]. Additionally, and last decade (Figure 3A). This phenological advance is pre- coarse talus slopes holding pockets of deeper soils in commu- dicted to be countered by an expected 13-day delay in flowering nities 3 and 4 might retain more water than is possible in the lower generated by the 11.6-cm reduction in annual precipitation be- two communities, and the highly fractured bedrock could provide tween the decades (Figure 3B; Table S1), yielding a predicted water storage in community 5 if refreshed by precipitation [61]. net advance of 28 days, only 8 days more than the observed These differences in physical features and evapotranspiration advance. The predicted magnitude of the effects of temperature rates (‘‘environmental mechanisms’’ sensu [59]) could explain Current Biology 30, 432–441, February 3, 2020 437 Figure 4. Change in Flowering Synchrony of Subpopulations Change in the interval between mean flowering times for subpopulations within communities 1 versus 2, 2 versus 3, 3 versus 4, and 4 versus 5 during the first decade (1984–1993) and last decade (2007–2016) of the survey period. Positive values indicate longer intervals (decreased synchrony), and negative values indicate shorter intervals (increased synchrony). The difference in flowering time for subpopulations on elevation band 2 versus 3 was significantly larger in the last decade than in the first, whereas the difference in flowering time for subpopulations on elevation band 4 versus 5 was significantly smaller in the last decade than in the first. Highest posterior density intervals (95%) are shown. See also Table S2. changing climatic cues. Our study pro- vides a novel view of how the timing of flowering is changing in a semi-arid ecosystem, an ecosystem type that is ex- pected to expand with continued climate why flowering shifted at different rates across the gradient, such change [63]. Ephemeral and intermittent stream communities as communities 1 and 2 shifting earlier at rates 3–7 times faster in semi-arid ecosystems contain high biodiversity and provide than did communities 3 and 4. In addition, the responses of com- the same ecological services as do true riparian areas [64]; munities 3 and 4 were shaped by significant interactions between these systems, and the biodiversity they contain, are highly temperature and precipitation (Table S1), and flowering in those threatened by climate change [65]. The phenological changes communities advanced more when conditions were both warmer driven by increasing temperature and decreasing precipitation and wetter. Though correlational, this result suggests that these in the xeroriparian habitat in our study system could be indica- interacting drivers can give rise to divergent phenologies in neigh- tive of what we can expect elsewhere in the United States, and boring communities. particularly in the Southwest. This study shows that ecosystem Given the high variability of precipitation and recurring drought responses to climate change will be both variable and complex, in the Southwest [34], it is not surprising that we did not detect a particularly in highly heterogeneous systems characterized by long-term linear trend in total annual precipitation (Figure S2B). high interannual climatic variability, highly variable topography, However, since the late 1990s, the study area has been impacted and high biodiversity. Short-term studies of only a few by long-term warm hydraulic drought that is characterized by not species might not show the extent of change occurring, espe- only precipitation deficits but also rising temperatures [33,34]. cially when baseline data are lacking, and could in fact produce The considerable increase in precipitation variability in the South- erroneous conclusions. Our findings demonstrate that commu- west could exacerbate the effects of drought by reducing growth nities and subpopulations occupying different microclimates and increasing mortality [35]. In addition, the increasing intensity are exhibiting remarkably different responses to changing of monsoon storms [36] will likely result in greater run-off, which climatic conditions. The differences in both magnitude and di- means less moisture available to vegetation. In the study area, rection of responses highlight how climate change will result several dominant species have declined in abundance, particu- in community reshuffling in the temporal dimension. larly Carnegiea gigantea (saguaro), Juniperus deppeana var. dep- peana (alligator juniper), Parkinsonia microphylla (foothills palo STAR+METHODS verde), Pinus discolor (border pinyon), Pinus ponderosa subsp. brachyptera (ponderosa pine), and Quercus arizonica (Arizona Detailed methods are provided in the online version of this paper white oak), as have nearly all annuals, particularly when cool-sea- and include the following: son rains are poor [62]. The result of these declines is likely d KEY RESOURCES TABLE increased soil temperatures due to an increase in bare ground d LEAD CONTACT AND MATERIALS AVAILABILITY and decreased cover. These drought conditions could explain d EXPERIMENTAL MODEL AND SUBJECT DETAILS some of the temporally nonlinear responses we detected, leading d METHOD DETAILS communities 3 and 4 to exhibit delayed flowering in the last B Phenological Dataset decade of the time series in relation to the first. B Climate Data By virtue of long-term, taxonomically extensive, and highly d QUANTIFICATION AND STATISTICAL ANALYSIS temporally resolved data that span a spatial gradient, we B Circular Statistics were able to detect divergent metacommunity-level and sub- B Metacommunity Shifts in Flowering Phenology population-level shifts in flowering phenology driven by 438 Current Biology 30, 432–441, February 3, 2020 B Climate Models 11. Rafferty, N.E., Bertelsen, C.D., and Bronstein, J.L. (2016). Later flowering is associated with a compressed flowering season and reduced reproduc- B Subpopulation Changes in Flowering Synchrony tive output in an early season floral resource. Oikos 125, 821–828. d DATA AND CODE AVAILABILITY 12. Inouye, D.W. (2008). Effects of climate change on phenology, frost dam- age, and floral abundance of montane wildflowers. Ecology 89, 353–362. SUPPLEMENTAL INFORMATION 13. Iler, A.M., Compagnoni, A., Inouye, D.W., Williams, J.L., CaraDonna, P.J., Anderson, A., and Miller, T.E. (2019). Reproductive losses due to climate Supplemental Information can be found online at https://doi.org/10.1016/j. change-induced earlier flowering are not the primary threat to plant popu- cub.2019.11.071. lation viability in a perennial herb. J. Ecol. 107, 1931–1943. 14. Pearse, W.D., Davis, C.C., Inouye, D.W., Primack, R.B., and Davies, T.J. ACKNOWLEDGMENTS (2017). A statistical estimator for determining the limits of contemporary and historic phenology. Nat. Ecol. Evol. 1, 1876–1882. N.E.R. acknowledges UCR initial complement funding and thanks Paul Nabity for insightful discussions and David Rafferty for assistance with coding to 15. Kudo, G., and Cooper, E.J. (2019). When spring ephemerals fail to meet extract mean differences in flowering time within species. C.D.B. acknowl- pollinators: mechanism of phenological mismatch and its impact on plant edges the University of Arizona Herbarium staff who identified or verified nearly reproduction. Proc. Biol. Sci. 286, 20190573. all the taxa in the flora. We thank Joel Sachs and four anonymous reviewers for 16. Parsche, S., Frund, J., and Tscharntke, T. (2011). Experimental environ- comments that greatly improved the manuscript. mental change and mutualistic vs. antagonistic plant flower-visitor interac- tions. Perspect. Plant Ecol. Evol. Syst. 13, 27–35. AUTHOR CONTRIBUTIONS 17. Rafferty, N.E., and Ives, A.R. (2011). Effects of experimental shifts in flow- ering phenology on plant-pollinator interactions. Ecol. Lett. 14, 69–74. C.D.B. collected the data, N.E.R. analyzed the data with input from J.M.D., 18. Gezon, Z.J., Inouye, D.W., and Irwin, R.E. (2016). Phenological change in a N.E.R. wrote the first draft of the manuscript, and all authors devised the study spring ephemeral: implications for pollination and plant reproduction. questions and edited the manuscript. Glob. Change Biol. 22, 1779–1793. 19. Rafferty, N.E., and Ives, A.R. (2012). Pollinator effectiveness varies with DECLARATION OF INTERESTS experimental shifts in flowering time. Ecology 93, 803–814. 20. Forrest, J.R. (2015). Plant-pollinator interactions and phenological The authors declare no competing interests. change: what can we learn about climate impacts from experiments and observations? Oikos 124, 4–13. Received: August 25, 2019 Revised: October 20, 2019 21. Burkle, L.A., Marlin, J.C., and Knight, T.M. (2013). Plant-pollinator interac- Accepted: November 25, 2019 tions over 120 years: loss of species, co-occurrence, and function. Published: January 2, 2020 Science 339, 1611–1615. 22. Bartomeus, I., Ascher, J.S., Wagner, D., Danforth, B.N., Colla, S., REFERENCES Kornbluth, S., and Winfree, R. (2011). Climate-associated phenological advances in bee pollinators and bee-pollinated plants. Proc. Natl. Acad. 1. Bradley, N.L., Leopold, A.C., Ross, J., and Huffaker, W. (1999). Sci. USA 108, 20645–20649. Phenological changes reflect climate change in Wisconsin. Proc. Natl. 23. Thomson, J.D. (2019). Progressive deterioration of pollination service de- Acad. Sci. USA 96, 9701–9704. tected in a 17-year study vanishes in a 26-year study. New Phytol. 224, 2. Fitter, A.H., and Fitter, R.S. (2002). Rapid changes in flowering time in 1151–1159. British plants. Science 296, 1689–1691. 24. Yang, L.H., and Rudolf, V.H. (2010). Phenology, ontogeny and the effects 3. Crimmins, T.M., Crimmins, M.A., and Bertelsen, C.D. (2010). Complex re- of climate change on the timing of species interactions. Ecol. Lett. 13, sponses to climate drivers in onset of spring flowering across a semi-arid 1–10. elevation gradient. J. Ecol. 98, 1042–1051. 25. Aldridge, G., Inouye, D.W., Forrest, J.R., Barr, W.A., and Miller-Rushing, 4. Crimmins, T.M., Crimmins, M.A., and Bertelsen, C.D. (2011). Onset of A.J. (2011). Emergence of a mid-season period of low floral resources in summer flowering in a ‘Sky Island’ is driven by monsoon moisture. New a montane meadow ecosystem associated with climate change. J. Ecol. Phytol. 191, 468–479. 99, 905–913. 5. Calinger, K.M., Queenborough, S., and Curtis, P.S. (2013). Herbarium 26. Theobald, E.J., Breckheimer, I., and HilleRisLambers, J. (2017). Climate specimens reveal the footprint of climate change on flowering trends drives phenological reassembly of a mountain wildflower meadow com- across north-central North America. Ecol. Lett. 16, 1037–1044. munity. Ecology 98, 2799–2812. 6. Diez, J.M., Iba´ n˜ ez, I., Miller-Rushing, A.J., Mazer, S.J., Crimmins, T.M., 27. Ison, J.L., Wagenius, S., Reitz, D., and Ashley, M.V. (2014). Mating be- Crimmins, M.A., Bertelsen, C.D., and Inouye, D.W. (2012). Forecasting tween Echinacea angustifolia (Asteraceae) individuals increases with their phenology: from species variability to community patterns. Ecol. Lett. flowering synchrony and spatial proximity. Am. J. Bot. 101, 180–189. 15, 545–553. 28. Weis, A.E., Nardone, E., and Fox, G.A. (2014). The strength of assortative 7. Willis, C.G., Ruhfel, B., Primack, R.B., Miller-Rushing, A.J., and Davis, C.C. mating for flowering date and its basis in individual variation in flowering (2008). Phylogenetic patterns of species loss in Thoreau’s woods are schedule. J. Evol. Biol. 27, 2138–2151. driven by climate change. Proc. Natl. Acad. Sci. USA 105, 17029–17033. 29. Waser, N., and Price, M. (1991). Outcrossing distance effects in 8. Rafferty, N.E., and Nabity, P.D. (2017). A global test for phylogenetic signal Delphinium nelsonii: pollen loads, pollen tubes, and seed set. Ecology in shifts in flowering time under climate change. J. Ecol. 105, 627–633. 72, 171–179. 9. Miller-Rushing, A.J., Høye, T.T., Inouye, D.W., and Post, E. (2010). The ef- 30. Morellato, L., Alberton, B., Alvarado, S.T., Borges, B., Buisson, E., fects of phenological mismatches on demography. Philos. Trans. R. Soc. Camargo, M.G., Cancian, L.F., Carstensen, D.W., Escobar, D.F., Leite, Lond. B Biol. Sci. 365, 3177–3186. P.T., et al. (2016). Linking plant phenology to conservation biology. Biol. Conserv. 195, 60–72. 10. Cleland, E.E., Allen, J.M., Crimmins, T.M., Dunne, J.A., Pau, S., Travers, S.E., Zavaleta, E.S., and Wolkovich, E.M. (2012). Phenological tracking en- 31. Halbritter, A.H., Alexander, J.M., Edwards, P.J., and Billeter, R. (2013). ables positive species responses to climate change. Ecology 93, 1765– How comparable are species distributions along elevational and latitudinal 1771. climate gradients? Glob. Ecol. Biogeogr. 22, 1228–1237. Current Biology 30, 432–441, February 3, 2020 439 32. Ko¨ rner, C. (2007). The use of ‘altitude’ in ecological research. Trends Ecol. 49. Franks, S.J., Weber, J.J., and Aitken, S.N. (2014). Evolutionary and plastic Evol. 22, 569–574. responses to climate change in terrestrial plant populations. Evol. Appl. 7, 123–139. 33. Arizona State Climate Office (2019). Drought. Global Institute of 50. Merila, € J., and Hendry, A.P. (2014). Climate change, adaptation, and Sustainability, Arizona State University, Tempe, AZ. https://azclimate. phenotypic plasticity: the problem and the evidence. Evol. Appl. 7, 1–14. asu.edu/drought/. 51. Miller-Rushing, A.J., Inouye, D.W., and Primack, R.B. (2008). How well 34. Woodhouse, C.A., Meko, D.M., MacDonald, G.M., Stahle, D.W., and do first flowering dates measure plant responses to climate change? Cook, E.R. (2010). A 1,200-year perspective of 21st century drought in The effects of population size and sampling frequency. J. Ecol. 96, southwestern North America. Proc. Natl. Acad. Sci. USA 107, 21283– 1289–1296. 52. Elzinga, J.A., Atlan, A., Biere, A., Gigord, L., Weis, A.E., and Bernasconi, G. 35. Dannenberg, M.P., Wise, E.K., and Smith, W.K. (2019). Reduced tree (2007). Time after time: flowering phenology and biotic interactions. growth in the semiarid United States due to asymmetric responses to Trends Ecol. Evol. 22, 432–439. intensifying precipitation extremes. Sci. Adv. 5, w0667. 53. Kottek, M., Grieser, J., Beck, C., Rudolf, B., and Rubel, F. (2006). World 36. Demaria, E.M.C., Hazenberg, P., Scott, R.L., Meles, M.B., Nichols, M., and map of the Ko¨ ppen-Geiger climate classification updated. Meteorol. Z. Goodrich, D. (2019). Intensification of North American monsoon rainfall as (Berl.) 15, 259–263. observed from a long-term high-density gauge network. Geophys. Res. 54. Wolkovich, E.M., Cook, B.I., and Davies, T.J. (2014). Progress towards an Lett. 46, 6839–6847. interdisciplinary science of plant phenology: building predictions across 37. Crimmins, T.M., Crimmins, M.A., Bertelsen, D., and Balmat, J. (2008). space, time and species diversity. New Phytol. 201, 1156–1162. Relationships between alpha diversity of plant species in bloom and cli- 55. Morellato, L.P.C., Alberti, L.F., and Hudson, I.L. (2010). Applications of cir- matic variables across an elevation gradient. Int. J. Biometeorol. 52, cular statistics in plant phenology: a case studies approach. In 353–366. Phenological Research, L.L. Hudson, and M.R. Keatley, eds. (Dordrecht: 38. Crimmins, T.M., Crimmins, M., and Bertelsen, C.D. (2013). Temporal pat- Springer), pp. 339–359. terns in species flowering in Sky Islands of the Sonoran Desert ecoregion. 56. Stein, A., Gerstner, K., and Kreft, H. (2014). Environmental heterogeneity In Merging Science and Management in a Rapidly Changing World: as a universal driver of species richness across taxa, biomes and spatial Biodiversity and Management of the Madrean Archipelago III, G.J. scales. Ecol. Lett. 17, 866–880. Gottfried, P.F. Ffolliott, B.S. Gebow, and L.G. Eskew, eds. (U.S. 57. Sundqvist, M.K., Sanders, N.J., and Wardle, D.A. (2013). Community Department of Agriculture, Forest Service, Rocky Mountain Research and ecosystem responses to elevational gradients: processes, mecha- Station), pp. 33–39. nisms, and insights for global change. Annu. Rev. Ecol. Evol. Syst. 44, 39. Crimmins, T.M., Crimmins, M.A., and Bertelsen, C.D. (2013). Spring and 261–280. summer patterns in flowering onset, duration, and constancy across a wa- 58. Wolkovich, E.M., Cook, B.I., Allen, J.M., Crimmins, T.M., Betancourt, J.L., ter-limited gradient. Am. J. Bot. 100, 1137–1147. Travers, S.E., Pau, S., Regetz, J., Davies, T.J., Kraft, N.J., et al. (2012). 40. Crimmins, T.M., Bertelsen, D.C., and Crimmins, M.A. (2014). Within-sea- Warming experiments underpredict plant phenological responses to son flowering interruptions are common in the water-limited Sky Islands. climate change. Nature 485, 494–497. Int. J. Biometeorol. 58, 419–426. 59. Chmura, H.E., Kharouba, H.M., Ashander, J., Ehlman, S.M., Rivest, E.B., 41. Iler, A.M., Høye, T.T., Inouye, D.W., and Schmidt, N.M. (2013). Long-term and Yang, L.H. (2019). The mechanisms of phenology: the patterns and trends mask variation in the direction and magnitude of short-term pheno- processes of phenological shifts. Ecol. Monogr. 89, e01337–e22. logical shifts. Am. J. Bot. 100, 1398–1406. 60. Whittaker, R., Buol, S., Niering, W., and Havens, Y. (1968). A soil and vege- 42. Iler, A.M., Høye, T.T., Inouye, D.W., and Schmidt, N.M. (2013). Nonlinear tation pattern in the Santa Catalina Mountains, Arizona. Soil Sci. 105, flowering responses to climate: are species approaching their limits of 440–450. phenological change? Philos. Trans. R. Soc. Lond. B Biol. Sci. 368, 61. Schwinning, S. (2010). The ecohydrology of roots in rocks. Ecohydrology 3, 238–245. 43. de Keyzer, C.W., Rafferty, N.E., Inouye, D.W., and Thomson, J.D. (2017). 62. Bertelsen, C.D. (2018). Thirty-seven years on a mountain trail: vascular Confounding effects of spatial variation on shifts in phenology. Glob. flora and flowering phenology of the Finger Rock Canyon watershed, Change Biol. 23, 1783–1791. Santa Catalina Mountains, Arizona. Desert Plants 34, 3–268. 44. Hendry, A.P., and Day, T. (2005). Population structure attributable to 63. Rajaud, A., and de Noblet-Ducoudre,  N. (2017). Tropical semi-arid regions reproductive time: isolation by time and adaptation by time. Mol. Ecol. expanding over temperate latitudes under climate change. Clim. Change 14, 901–916. 144, 703–719. 45. Vitasse, Y., Signarbieux, C., and Fu, Y.H. (2018). Global warming leads to 64. Levick, L., Fonseca, J., Goodrich, D., Hernandez, M., Semmens, D., more uniform spring phenology across elevations. Proc. Natl. Acad. Sci. Stromberg, J., Leidy, R., Scianni, M., Guertin, D.P., Tluczek, M., et al. USA 115, 1004–1008. (2008). The ecological and hydrological significance of ephemeral and intermittent streams in the arid and semi-arid American Southwest. 46. Prevey, J., Vellend, M., Ru¨ ger, N., Hollister, R.D., Bjorkman, A.D., Myers- https://www.epa.gov/sites/production/files/2015-03/documents/ephemeral_ Smith, I.H., Elmendorf, S.C., Clark, K., Cooper, E.J., Elberling, B., et al. streams_report_final_508-kepner.pdf. (2017). Greater temperature sensitivity of plant phenology at colder sites: implications for convergence across northern latitudes. Glob. Change 65. Perry, L.G., Andersen, D.C., Reynolds, L.V., Nelson, S.M., and Shafroth, P.B. (2012). Vulnerability of riparian ecosystems to elevated CO2 and Biol. 23, 2660–2671. climate change in arid and semiarid western North America. Glob. 47. Wright, K.W., Vanderbilt, K.L., Inouye, D.W., Bertelsen, C.D., and Change Biol. 18, 821–842. Crimmins, T.M. (2015). Turnover and reliability of flower communities in 66. Brown, D.E. (1994). Biotic communities: Southwestern United States and extreme environments: insights from long-term phenology data sets. Northwestern Mexico (Salt Lake City: University of Utah Press). J. Arid Environ. 115, 27–34. 67. PRISM Climate Group Oregon State University (2017). created 25 Oct. 48. Anderson, J.T., Inouye, D.W., McKinney, A.M., Colautti, R.I., and Mitchell- http://prism.oregonstate.edu. Olds, T. (2012). Phenotypic plasticity and adaptive evolution contribute to advancing flowering phenology in response to climate change. Proc. Biol. 68. Daly, C., Halbleib, M., Smith, J.I., Gibson, W.P., Doggett, M.K., Taylor, Sci. 279, 3843–3852. G.H., Curtis, J., and Pasteris, P.P. (2008). Physiographically sensitive 440 Current Biology 30, 432–441, February 3, 2020 mapping of climatological temperature and precipitation across the 72. Cremers, J., Mulder, K.T., and Klugkist, I. (2018). Circular interpretation of conterminous United States. Int. J. Climatol. 28, 2031–2064. regression coefficients. Br. J. Math. Stat. Psychol. 71, 75–95. 69. Agostinelli, C., and Lund, U. (2017). R package ‘circular’: circular statistics 73. R Core Team (2019). R: A language and environment for statistical (version 0.4-93). URL. https://r-forge.r-project.org/projects/circular/. computing. R Foundation for Statistical Computing, Vienna, Austria. URL. https://www.R-project.org/. 70. Cremers, J. (2018). bpnreg: Bayesian projected normal regression models for circular data. R package version 1.0.0. URL. https://CRAN.R-project. 74. Pewsey, A., Neuhauser, M., and Ruxton, G.D. (2013). Circular Statistics in org/package=bpnreg. R (New York: Oxford University Press). 71. Cremers, J., and Klugkist, I. (2018). One direction? A tutorial for circular 75. Iler, A.M., Inouye, D.W., Schmidt, N.M., and Høye, T.T. (2017). Detrending data analysis using R with examples in cognitive psychology. Front. phenological time series improves climate-phenology analyses and re- Psychol. 9, 2040. veals evidence of plasticity. Ecology 98, 647–655. Current Biology 30, 432–441, February 3, 2020 441 STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Deposited Data Flowering phenology data This paper; Mendeley data Mendeley data (https://doi.org/10.17632/k6p34z78x9.1) Climate data This paper; Mendeley data Mendeley data (https://doi.org/10.17632/k6p34z78x9.1) LEAD CONTACT AND MATERIALS AVAILABILITY Further information and requests should be directed to and will be fulfilled by the Lead Contact, Nicole Rafferty (nicole.rafferty@ucr. edu). The datasets used in this study have been deposited to Mendeley data (https://doi.org/10.17632/k6p34z78x9.1). This study did not generate any new reagents. EXPERIMENTAL MODEL AND SUBJECT DETAILS The data come from systematic surveys by one coauthor (CDB) of all plant species and infraspecific taxa (hereafter ‘‘species’’) in flower along a trail in Finger Rock Canyon ascending to Mt. Kimball in the Santa Catalina Mountains of Arizona, USA (Figure S3). Although the canyon represents less than 1% of the area of the Santa Catalina range, 45% of the known plant taxa in the range have been found there [62]. In 8.05 km, the trail ascends from 945-2212 m, which was partitioned at the beginning of data collection by CDB into five elevation bands that include six vegetative associations in five biotic communities (based on [66]): 1) 945-1079 m (desert scrub, riparian scrub); 2) 1079-1372 m (desert scrub, scrub grassland); 3) 1372-1671 m (scrub grassland, oak woodland); 4) 1671-1939 m (oak-pine woodland); 5) 1939-2212 m (oak-pine woodland, pine forest; Figure S3). METHOD DETAILS Phenological Dataset Every species seen in anthesis (angiosperms) or releasing pollen (gymnosperms), together referred to as ‘‘flowering,’’ was recorded for each community along each 1.6-km-long trail segment on every survey. During the first nine years of data collection, a period characterized by above-normal precipitation, data were collected an average of 30 days per year, with at least two surveys per month during the growing seasons. Subsequently, data were collected an average of 50 days per year, nearly weekly. Because our analyses use all records of flowering and focus on mean flowering dates, this change in sampling frequency should not bias our estimates. Surveys were completed throughout the year with approximately 8% of the total number of surveys being completed each month of the year. The 33-year survey period (1984-2016) considered here comprises 169,030 observations collected during 1,639 surveys. Of the 590 species, 117 were observed in only one community, 140 were observed in two communities, 168 were observed in three communities, 100 were observed in four communities, and 65 were observed in all five communities. Additional details about the data collection protocol and transect can be found in Crimmins et al. [39] and Bertelsen [62]. In particular, Bertelsen [62] gives for each species the years flowering and the number of flowering observations per elevation band and month. Climate Data The primary source of climate data was the Parameter-elevation Regressions on Independent Slopes Model (PRISM) database [67], supplemented by on-site rain gauges. Gauges (Tru-Chek) were installed by one of us (CDB) in 2007 to obtain data specific to three locations: at 957 m (near the base of the transect), 1459 m (approximately midway up the transect), and 2206 m (near the peak of the transect). Each gauge was checked on average four times per month during 2007-2012 and 2014-2016, and mineral oil was used to prevent evaporation. Temperature patterns and precipitation data for years the gauge data were not available were extracted from 4 km PRISM cells that include the GPS coordinates of the gauges. PRISM data incorporates factors such as location, elevation, and topography in a climate-elevation regression for each grid cell [68]. Although two of the gauges are located within the same PRISM cell, GPS coordinates within the cell produce different values based on elevation. Monthly PRISM cell and rain gauge precipitation records are highly correlated for each of the three locations (r = 0.85-0.89). Thus, monthly temperature and precipitation data were extracted for the same approximate elevations as the gauge locations to create three elevation-specific climate predictors of flower- ing phenology at low-, mid-, and high-elevation. Based on detailed knowledge of the aspect and topography of each elevation band, long-term observation of weather patterns, and vegetation responses to short-term climatic events [see 62], we used low-elevation temperature and precipitation data to represent climate variables for communities 1 and 2, mid-elevation data to represent climate variables for community 3, and high-elevation data to represent climate variables for communities 4 and 5. Briefly, communities 1 and 2 have significant southern exposure, likely resulting in higher evapotranspiration, particularly since precipitation is less than at higher e1 Current Biology 30, 432–441.e1–e3, February 3, 2020 elevations. Community 3 is situated in the deepest and narrowest portion of the canyon, measured from the top of the ridges forming the canyon to the bottom of the drainage; temperatures are likely moderated by cold air drainage and the largely southwestern expo- sure. Community 4 has more continuous cover and a largely northwestern exposure, whereas community 5 has considerable cover but with extensive areas of exposed bedrock. The higher precipitation received in communities 3-5, the amount of cover, and expo- sures would likely lessen evapotranspiration. We regressed mean daily temperature and total annual precipitation against year to test whether these climate variables have changed over the survey period (1984-2016) and since 1930, when a sufficient number of nearby weather stations were available to provide reliable data for the study area [35]. In the initial regression models, we included the interaction between year and elevation to test whether slopes differed by elevation, in which case we fit separate regressions per elevation. If the interaction term was not significant, we fit a linear mixed-effects model with elevation as a random effect to allow intercepts to vary by elevation and compared the fit of models with and without this random effect (conditional versus marginal R ). QUANTIFICATION AND STATISTICAL ANALYSIS Circular Statistics Data from thirty years were used in the statistical analysis; 2004, 2005, and 2013 were excluded because surveys occurred irregularly during those years. We converted all survey dates to day of year (doy). For all circular statistics, we converted doy to radians and used the R packages ‘‘circular’’ version 0.4-93 [69] and, to construct circular mixed-effects models, ‘‘bpnreg’’ version 1.0.0 [70]. Additional details regarding circular mixed-effects models, their formulation, and interpretation can be found in Cremers and Klugkist [71] and Cremers et al. [72]. Briefly, for the circular mixed-effects models, statistical significance of continuous predictors was gauged by whether the 95% highest posterior density (HPD) lower and upper bounds for circular model coefficients included zero (not signif- icant); significance of categorical variables was judged by whether the 95% HPD intervals for both component I and component II linear coefficients included zero (not significant; [72]). For circular differences between variables, significance was determined by the proportion of iterations that were negative. Models were compared using the deviance information criterion (DIC and DIC ) alt and the Watanabe-Akaike information criterion (WAIC and WAIC ), which both reward better-fitting models while penalizing model 1 2 complexity [70]. All analyses were conducted in R version 3.5.2 [73]. Metacommunity Shifts in Flowering Phenology To examine how flowering phenology at the metacommunity level has changed over time while holding elevation constant, we compared all dates on which flowers of any species were observed within a given community from 1984-2016. This analysis allowed us to determine whether communities within each elevation band exhibited trends of advanced or delayed flowering across all years. We fitted additional models that included a quadratic term for year to test for nonlinear changes in flowering phenology over the full time series. We were also interested in examining trends in the early versus later years of the survey period with the expectation that any phenological changes in response to climate change would be most apparent when comparing the two end segments of the time series. In particular, phenological effects of the drought that began in the late 1990s are likely to be apparent in the latest years of the time series. Therefore, we also examined how flowering time changed per community between the first and last decades of the data- set (i.e., 1984-1993 versus 2007-2016), which additionally enabled us to visualize any shifts with circular plots. Because these ana- lyses use all available flowering records for each community, the various flowering distributions (including any bimodal distributions) are aggregated when mean flowering times are calculated. For each community, we constructed circular mixed-effects models with doy of flowering (in radians) as the response and year (continuous), year + year , or decade (categorical) as the predictor(s), with spe- cies identity included as a random effect to account for repeated observations of the same species over time. Year was centered at zero to aid interpretation of model coefficients. We also verified that communities exhibited elevation-band-specific responses in phenology by fitting a circular mixed-effects model with year, elevation band (continuous), and the interaction between year and elevation band as predictors (and species identity as a random effect). The bpnme () function within the bpnreg package uses a Bayesian approach and Markov chain Monte Carlo (MCMC) samplers to estimate model parameters [71]. For each model, we ran 10,000 iterations with a burn-in period of 1,000 iterations and no lag because there was minimal auto-correlation detected in the MCMC chains. Because the model with the interaction between year and elevation band as a predictor was very computationally intensive, we ran only 1,000 iterations with a burn-in period of 100 iterations for that model. We inspected traceplots to verify that models had converged, which is the method recommended by the package authors [70]. For models with continuous predictors, we report the ‘slope at the mean’ (SAM) circular model coefficients because they are the least biased [72]. To test whether the con- centration of flowering times (i.e., the spread of flowering times throughout the calendar year) changed for a given community be- tween decades, we used Wallraff’s test for a common concentration [74]. Climate Models To test whether climate variables were related to observed shifts in flowering time within each community, we used circular mixed- effects models with doy of flowering as the response and mean daily temperature and total precipitation during the 12-month window preceding and including the month in which flowering was observed as predictors, with species identity as a random effect. In our initial models, we also included the interaction between temperature and precipitation. These climate variables are specific to each flowering record and capture temperature and precipitation conditions for a one-year period before each observation, regardless of the calendar date of flowering. To test for temporal autocorrelation within each time series (precipitation, temperature, and flowering Current Biology 30, 432–441.e1–e3, February 3, 2020 e2 phenology), we used Ljung-Box tests with a lag of one year [75]. Temperature and precipitation variables were centered at zero to aid interpretation of model coefficients. As before, we ran 10,000 iterations per model with a burn-in period of 1,000 iterations and no lag because there was minimal auto-correlation detected in the MCMC chains. We inspected traceplots to verify that MCMC chains had converged [70]. For these models, we report the SAM circular model coefficients [72]. Subpopulation Changes in Flowering Synchrony To analyze within-species shifts in flowering phenology across time and space, we limited the dataset to only those species that: (i) occurred in at least two adjacent communities along the transect (i.e., communities 1 and 2; communities 2 and 3; communities 3 and 4; communities 4 and 5), (ii) had been observed flowering during at least four years in the first decade (1984-1993) and four years in the last decade (2007-2016), and (iii) had been observed at least four times per community in each year flowering was documented. The resulting dataset comprised 128 unique species with subpopulations in two or more adjacent communities. We then calculated the mean doy of flowering (in radians) per species per community per decade and took the difference between these means for each community per decade (e.g., we subtracted the mean doy of flowering for a given species in community 2 in the first decade from the mean doy of flowering for the same species in community 1 in the first decade; we repeated this process for the last decade). These values indicate how much the mean flowering times of subpopulations differed in the first versus the last decade of the dataset and provide a way to test for changes in subpopulation-level synchrony. These differences were first converted back to a circular variable so that values corresponded to the number of days in radians on the unit circle, measuring counterclockwise from 0 radians if the differences were positive and clockwise from 2p radians if the differences were negative. We then used these differences as our response variable in circular mixed-effects models with decade as the predictor and species identity as a random effect. We used the same model specifications as in previous models. DATA AND CODE AVAILABILITY The phenology and climate datasets used in this study are available at Mendeley data (https://doi.org/10.17632/k6p34z78x9.1). e3 Current Biology 30, 432–441.e1–e3, February 3, 2020

Journal

Current BiologyUnpaywall

Published: Feb 1, 2020

There are no references for this article.