Phenology, the timing of biological events, from bird migrations to flowering, has become a ‘leading indicator’ of climate change impacts (Parmesan and Yohe 2003). Climate-induced shifts in phenology can negatively affect individual fitness if the new timing of a biological event does not correspond to when environmental conditions are optimal (Willis et al. 2008). Changes in the timing of growth and reproduction may ultimately be necessary to track the rapid pace of climate change (Cleland et al. 2012). Evaluations of how climate and phenology are coupled, and how they have changed through time, can provide insight on the potential for organisms to persist under climate change. Understanding whether species respond in similar or divergent ways to different aspects of climate is essential to understanding future community composition, trophic dynamics, and ecosystem function.

Decades of phenological research has demonstrated the practicality and utility of linking ground-based phenological and meteorological measurements. Many of these past studies have focused primarily on temperature and the accumulation of growing degree days as phenological drivers (e.g., Basler 2016). This research has demonstrated that warming temperatures and a more rapid accumulation of growing degree days (AGDD) have accelerated the phenology of many species. Growing degree days provide a measure of heat accumulation in the spring, defined as the number of degrees by which average daily temperatures exceed some baseline (e.g., freezing). They are accumulated (summed) daily following a predefined start date. AGDD has proven to be an effective predictor of phenological transitions in plant, insect, and other animal species (Cayton et al. 2015; Crimmins et al. 2017) leading to the development of gridded phenological indices that have been used to forecast and map the onset of spring (Schwartz 1997; Schwartz et al. 2013; Ault et al. 2015) and, more recently, species-level phenological forecasts (Taylor and White 2020). However, multiple aspects of plant biology, from phenology to growth and mortality, are also influenced by the balance between the timing and magnitude of precipitation and atmospheric demand for moisture (e.g., Novick et al. 2016; Choat et al. 2018). Understanding the windows of time in which precipitation, atmospheric demand for moisture, and temperature-related variables (e.g., AGDD) affect plant phenology is essential, and can reveal how changes in moisture availability mediated by climate change may influence ecosystems.

Although temperature-based approaches to phenology have been successful, a rich history of phenological research has shown that rainfall and moisture availability also play an important role in the timing of phenological events from tropical rainforests to the water-limited Sonoran Desert of North America (Gentry 1974; Reich and Borchert 1984; Ashton et al. 1988; Fox 1990; Bowers and Dimmitt 1994; Lasky et al. 2015). Nearly three decades ago, Bowers and Dimmitt (1994) showed that the flowering of six dominant woody species is triggered by rain and photoperiod, after sufficient accumulations of growing degree days. Since that time, the role of rainfall has received increasing attention in phenological research in dryland ecosystems (Mazer et al. 2015; Park and Mazer 2018; Renzi et al. 2019; Elmendorf et al. 2019). Another important aspect of moisture availability is the “thirst” of the atmosphere, or atmospheric demand. The role of atmospheric demand is only now beginning to appear in the phenological literature, but recent research has shown that increases in atmospheric demand can alter, often delaying, plant phenology (Adams et al. 2015; Wion et al. 2020).

Many of the climatic factors that influence plant performance are rapidly changing in the southwestern U.S. Temperatures across the Sonoran Desert have risen by 1–3°C in the last century (Munson et al. 2012) and the southwestern U.S. is expected to rise 3–6°C by 2100 (IPCC 2013). This warming has extended the length of the frost-free season (Osland et al. 2021) and – in conjunction with lower relative humidity – has been associated with prolonged drought conditions since 2000 (Weiss and Overpeck 2005; Weiss et al. 2009). Indeed, in the southwestern U.S., temperature and vapor pressure deficits (VPD) are increasing faster than in other regions (Driscoll et al. 2020). In the Mojave Desert, for example, there have been substantial increases in both mean annual temperature and mean daily maximum VPD, alongside decreases in total annual precipitation, resulting in much lower ecologically available water (Driscoll et al. 2020). In general, there is broad consensus that aridity will continue to increase in the Desert Southwest due to both decreasing precipitation and increasing atmospheric demand as we move into the mid to late 21st century (Seager and Vecchi 2010; Cayan et al. 2010; Udall and Overpeck 2017).

How will phenology in the Desert Southwest respond to recent and ongoing climate change? Bowers and Dimmitt (1994) showed that six dominant woody species in the Sonoran Desert have a spectrum of phenological strategies that vary in their relation to photoperiod, temperature, and rainfall. These findings are part of a growing body of research demonstrating variation in phenological responses to climate change across species (Fitter and Fitter 2002; Parmesan 2007; Sherry et al. 2007; Crimmins et al. 2010; CaraDonna et al. 2014; CaraDonna and Inouye 2015). Here, we revisit the work of Bowers and Dimmitt (1994), focusing on many of the same species, but with a more complete phenology record and new phenological modeling techniques. With this approach, our objectives were to determine (1) how spring flowering time has changed over the last three decades, (2) which aspects of climate have influenced the probability and timing of flowering, and (3) how climatic factors that have influenced flowering phenology historically are forecasted to change into the future in the Sonoran Desert. We focus our attention on spring flowering time because of the potential for phenological advances due to warming, phenological triggers due to moisture availability from cool-season frontal storms, and the higher consistency of floral output and pollinator activity across species at this time of year. The climatic factors we focused on were selected because they are biologically meaningful drivers of flowering phenology in the Sonoran Desert, and include temperature variables (growing degree days, vapor pressure deficit, and number of freezes) and precipitation variables at multiple time scales (2-week rainfall, number of days since rainfall, antecedent seasonal rainfall, and the Pacific Decadal Oscillation Index). We contextualize our phenological modeling results with existing climate projections to understand how the patterns we observed might play out to the end of the century.

Methods

Study Species

We examined the flowering phenology of six dominant drought-deciduous woody perennials in the Arizona Upland subdivision of the Sonoran Desert for which we had an abundance of phenological data (Table 1). Drought-deciduous shrubs and trees represent an important functional group within warm deserts (Shreve 1951). Although as a group they have adopted the drought-evading strategy of shedding leaves in response to increasing water stress, the morphological and physiological diversity within this functional group is immense and has resulted in a large diversity of phenological strategies.

Table 1.

Phenological Windows and Functional Traits of Focal Species. The start date (mm/dd) marks the beginning of the window within which we might expect to see flowering. The end date (mm/dd) is used to ensure that monsoon season events, which occur after phenological events during the usual cool-season period, are not accidentally included in the model. Our focus throughout is on “cool” (winter/spring) flowering events as opposed to monsoon season events. For context, we also provide additional growth habit and life history information.

Two of our species, Parkinsonia microphylla (Fabaceae) and Olneya tesota (Fabaceae) are trees, thought to be deep-rooted (Cannon 1911; Phillips 1963; Canadell et al. 1996), which flower between April–June. While O. tesota drops its leaves only once per year, just as flower buds appear, P. microphylla can remain leafless most of the year by relying on its photosynthetic trunk and branches (Smith 1997).

Vachellia constricta (Fabaceae), Whitethorn Acacia, is a drought-deciduous shrub, sometimes reaching the stature of a small tree, which blooms in late spring (May–June) and again in summer/autumn (July–October) (McGinnies 1983). It is the only one of our species which has evolved to take advantage of the bimodal distribution of rainfall in the Sonoran Desert by flowering in both the spring and summer.

Fouquieria splendens (Fouquieriaceae), Ocotillo, is a drought-deciduous, stem-succulent shrub (Killingbeck 2019), the only one of our species to have adopted both deciduousness and stem succulence in response to highly limited and highly variable water availability. In the Sonoran Desert, it blooms primarily in spring (March–May; Waser 1979), but massive autumn flowering can occur under some circumstances (Felger 1980), and winter (January–February) flowering is possible.

Encelia farinosa (Asteraceae), Brittlebush, and Ambrosia deltoidea (Asteraceae), Triangle-leaf Bursage, are drought-deciduous semi-woody subshrubs. Encelia farinosa blooms mainly in spring (February–May), and given enough rain and lack of frost, from October–January as well. Ambrosia deltoidea is one of the most abundant species in the Arizona Upland subdivision of the Sonoran Desert. It typically flowers from February through early April. This is the only wind-pollinated species within our focal taxa.

Of the original six species in Bowers and Dimmitt (1994) we excluded Larrea tridentata, Creosote Bush, which was found in flower in every month of the year during the 36 years of phenology observations collected at our study site. We sought to model the onset of cool (winter/spring) season flowering for each species, but defining a consistent cool season flowering phenology window for L. tridentata proved difficult and cast doubt on our ability to produce reliable results for this species.

Model and Covariates

Our models of phenological transitions included weather variables as drivers, which we calculated using Daymet, a 1 km resolution gridded data product that provides estimates of daily weather parameters, including minimum and maximum temperature, precipitation, day length, and the partial pressure of water vapor (Daymet V3; Thornton et al. 2016). Daymet data were obtained using Google Earth Engine (Gorelick et al. 2017). We developed precipitation-based covariates to consider various time scales including rolling 2-week rainfall sums, the number of days since the most recent rainfall event (≥5 mm; hereafter ‘rainless days’), and antecedent rainfall – i.e., total precipitation during the monsoon season (June 15–September 30) preceding the early, cool-season flowering that we focused on modeling in this study. We developed temperature-based covariates as accumulated growing degree days (over a 10°C baseline), and the number of freezes (the total number of days with minimum temperatures below 28°F (-2.22°C), a threshold commonly used to define “hard” freezes). We developed variables to represent vapor pressure deficit as a 2-week rolling mean. We derived VPD from Daymet’s partial pressure of water vapor (E_a) by first calculating saturated vapor pressure (E_s) with the appropriate Arden Buck equation (Buck 1981), using the formula

where T is the mean of minimum and maximum temperature in degrees Celsius. We then calculated VPD as E_s – E_a. Finally, we include variables to represent a large-scale climate oscillation that controls much of the decadal-scale variation in moisture in the region, the Pacific Decadal Oscillation (PDO) with the PDO Index (Zhang et al. 1997; Mantua et al. 1997).

We modeled daily probabilities of flowering and, by extension, the distribution of flowering times, for each of these species using covariates representing current and antecedent temperature, precipitation, and vapor pressure deficit, and the PDO index. All modeling was performed using the R package tempo (Landau and Zachmann 2019). The landing page for the software provides a mathematical description of the model. This Bayesian model predicts the occurrence of phenological state transitions (e.g., from non-flowering to flowering) while accommodating censoring, a circumstance in which the value of an observation is only partially known. Phenological data sets are often censored because the event of interest is known only to have occurred sometime between visits to a site, but the precise time of the event is unknown. This is often the case in phenological observations, where technicians may only be able to visit sampling sites intermittently, for example, once per week. The Bayesian model accounts for censoring by imputing event (i.e., flowering) times between bounds provided by – in our case, approximately weekly – visitation dates, and avoids the bias encountered when censoring is ignored.

Another important feature of the model is that it estimates phenological transitions on a daily basis, which allows us to represent temporally heterogeneous transition probabilities and their drivers. In other words, the chances an event is observed can increase when conditions are favorable, and subsequently decrease, at multiple times during the growing season. The shift in focus from modeling event time to a model of event occurrence has several important consequences (Clark et al. 2014; Diez et al. 2014; Elmendorf et al. 2019). One important strength of such an approach is that environmental variability is allowed to enter a model of event occurrence in an ecologically sensible way. For our data set, the complete daily history of rainfall events enters the model instead of seasonal or annual rainfall averages, which allows the probability of flowering to rise and fall over time. In contrast, models of event times require collapsing environmental variability into a single measure – a single summary measure per event, which ultimately results in less information entering model predictions. These models also do not allow for predictions of daily flowering probabilities, and by extension, they do not allow predictions of the probability that flowering will occur within a given time window. The model we used, however, does accommodate these sorts of predictions.

A bimodal precipitation regime has driven the evolution of two flowering seasons in the Sonoran Desert, with a cool (winter/spring) season and a warm (summer) season. We sought to model first flowering for each species during the cool (winter/spring) season. All time-to-event measures, as well as covariates, were developed using a phenological window that was unique to each species (Table 1). The phenological window for each species spanned the earliest to latest date at which onset of cool season flowering could conceivably occur. For instance, because cool season flowering for E. farinosa can occur as early as October, its phenological window spanned the transition between calendar years.

We ran a series of models, which all include AGDD and a quadratic term for AGDD, which allowed the probability of an event to increase and subsequently decrease within a growing season, as opposed to simply monotonically increasing. In addition to AGDD, we included rolling 2-week total precipitation (e.g., the value for a given day is the total precipitation that occurred in the preceding two weeks) and rolling 2-week mean vapor pressure deficit by default. We formed the rest of the models by including all possible combinations of the remaining covariates — number of freezes, rainless days, antecedent monsoonal precipitation, and the PDO index. All models were run with sufficient MCMC iterations to ensure convergence and appropriate characterization of posterior distributions. We used vague priors for all model parameters. We selected the “best” model for each species using the deviance information criterion (DIC; Spiegelhalter et al. 2002), excluding from consideration models that did not pass posterior predictive checks (Conn et al. 2018). We characterized trend by predicting mean flowering time for each year, followed by a regression of expected flowering times on year for each species.

We visualized the sensitivity to climate and the relative importance of each covariate in the “best” model in two ways. First, we shifted each covariate up or down to ± 30% of its observed range, while keeping all other covariates at their daily means, to measure the average departure of phenological transition time from the grand mean resulting from changes to covariates. Evaluating predictions conditional on covariates over some range (their native range or increases or decreases in their extremes) is a useful approach for understanding the relative importance of each covariate in a model. Second, we considered the two years representing the extreme ends of each covariate for each species in our historical record, with all other covariates held at their daily means, to visualize the effect on the distribution of flowering times.

Results

Climate conditions around the Arizona-Sonora Desert Museum are projected to change markedly by 2100 (Fig. 1). Gridded climate products show that mean temperature has been steadily rising in the historical record, and by 2100, is projected by the high-emissions RCP 8.5 scenario to be nearly 27°C, which is almost 6°C warmer than mean temperatures in the region in the 1980s (Fig. 1A). Annual rainfall has been declining historically, from ∼320 mm in the 1980s to ∼260 mm in the 2010s, however future projections range from either an increase of ∼280 mm or a decrease of ∼120 mm compared to the 1980s (Fig. 1B). The mean of future projections suggests a future average of 310 mm of annual rainfall, and thus little change from the 1980s. Like temperature, vapor pressure deficit has also risen in the historical period, and future projections show a clear rise in VPD, from just over 2 kPa in the 1980s to a range across all models of 2.5–3.5 kPa by 2100 (Fig. 1C). There is strong agreement across climate models for an increase of >1 kPa in VPD in this region (Ficklin and Novick 2017). These historical trends in annual summaries were reflected in the daily or weekly covariate data used in this study. For example, the actual AGDD data for each day of year for every year (1980–2019, Fig. 2), showed more rapidly accumulating growing degrees early in the year and overall higher accumulations by the end of the year in more recent years. Data at this resolution – days, weeks, or months – gave us an ability to model intraannual phenological forcings with much greater precision, while still reflecting the same overall trends in climate change visible in the annual summaries.

Fig. 1.

Annual summaries of historical and projected future climate conditions for the grid cell surrounding the Arizona-Sonora Desert Museum (32.2434°, -111.1672°). Gridded historical and future climate data were retrieved from the Climate Toolbox. Historical data are from the gridded surface meteorological data set, gridMET, which covers the contiguous U.S. at a 4 km spatial resolution from 1979–present. Projected future climate conditions include data from 20 climate models (the semi-opaque lines extending from 2020 to 2100) and one scenario (RCP 8.5) downscaled to a ∼4 km resolution across the U.S. for compatibility with the gridMET data using MACAv2-METDATA, version 2 (Abatzoglou and Brown 2012). The three panels show mean annual temperature (the average of maximum and minimum annual temperatures) (A), total annual precipitation (B), and mean vapor pressure deficit (C).

Fig. 2.

Accumulated Growing Degree Day (AGDD) data for each day of year for every year (1980–2019) at the Arizona-Sonora Desert Museum, computed using temperature information from Daymet. These daily measures provide an example of the type of climate/weather variable used in the phenological models developed for each species. The AGDD summaries seen here were computed from Jan 1 to Dec 31 of each year using a 10 °C baseline temperature.

Mean dates of flowering onset varied across the species in our study, from January 28th in E. farinosa, to May 9th in O. tesota (Table 2). Across our study period, five of six species appeared to exhibit an advance in the onset of flowering, ranging from 1.05 to 4.39 days earlier per decade (Table 2, Fig. 3,  Table S1 (473_madr-68-04-18_s01.docx)). One species, E. farinosa, showed an average predicted advance of 0.54 days per decade, but with high uncertainty. Our model predicted a ∼42% probability that E. farinosa was actually experiencing phenological delay over time (Table 2, Fig. 3,  Table S1 (473_madr-68-04-18_s01.docx)). The pattern of phenological advance seen among the rest of the species was more certain, with chances of phenological advance ranging from ∼84% for V. constricta to 100% for O. tesota and A. deltoidea (Table 2).

Table 2.

Summaries of Phenological Trends According to the Trend Lines in Fig. 3. Trends were estimated using a truncated normal linear regression of mean flowering onset dates on year. The expected mean flowering time in 1982, 2018, and over all years is reported to provide a general indication of month and day of first flower over the study period (conditional on the estimated trend). Also included are the shift, in days per decade (with 95% highest posterior density intervals provided parenthetically) and chances of a non-zero trend.

Fig. 3.

Median (points), 50% highest posterior density interval (HDI; inner thick error bars), and 95% HDI (outer thin bars) of predicted mean flowering times (y-axis) by species by year (x-axis). Note that the y-axes differ for each species because the phenological window for cool-season flowering for each species is different. Black error bars correspond to years during which a given species was observed. Gray bars denote unsampled years. The semi-opaque gray lines beneath individual mean flowering times provide an estimate of trend – i.e., phenological advance (decreasing slopes) or delay (increasing slope). Trend results were developed using a truncated normal linear regression of mean flowering times on year at each iteration of the MCMC algorithm used to fit the model for each species, which appropriately propagates uncertainty from the phenological model into the estimates of trend. The spread in gray lines represents uncertainty in the trend. For more complete information on trend and associated uncertainty, see Table 2.

Phenological models showed that the effect of AGDD, and its square (AGDD raised to the second power, which is done to model quadratic effects), were the most important factors for the onset of flowering across all species, with large coefficients (Fig. 4,  Fig. S1 (473_madr-68-04-18_s01.docx)). All species showed an advancing trend (earlier flowering) with increasing AGDD (Table 3, Fig. 4). Rolling 2-week rainfall advanced the flowering of P. microphylla and F. splendens and delayed the flowering of the remaining four keystone species (Table 3,  Fig. S1 (473_madr-68-04-18_s01.docx)). Rolling two-week mean VPD was important for every species, advancing flowering for O. tesota and E. farinosa and delaying the flowering for the other species ( Fig. S2 (473_madr-68-04-18_s01.docx)). Rainless days, antecedent monsoonal rainfall, the number of freezing events, and the PDO were only important for some species. Specifically, PDO advanced flowering in A. deltoidea and E. farinosa, and number of freezing events delayed flowering in O. tesota and F. splendens , while rainless days and antecedent monsoonal had mixed effects on flowering for the species they affected (O. tesota, A. deltoidea, and E. farinosa;  Fig. S2 (473_madr-68-04-18_s01.docx)). Divergent responses (advance vs. delay) across species to moisture-based climate variables was also apparent when we shifted the climate variables by ± 30% each day relative to that day's mean (Table 4,  Fig. S2 (473_madr-68-04-18_s01.docx)), and when we compared models representing the extreme ends of a particular climate variable while holding all other variables constant at their daily means (Fig. 5).

Fig. 4.

The expected influence of accumulated growing degree days (AGDD) on the distribution of flowering times in the two years representing the extreme ends of AGDD for each species. The lowest (dark grey lines) and highest (light grey lines) AGDD years varied by species because the phenological window within which AGDD was summarized varied by species, and the years for which we had data also varied by species. These results are based on the medians of all coefficient values reported in  Fig. S1 (473_madr-68-04-18_s01.docx). A with all other covariates held at their daily means. We show relative time to flowering (as opposed to date; mm/dd) within each species' phenological window along the x-axis so that effect size is more apparent. See Table 1 for detail on species' phenological windows. All species show an advancing trend (earlier flowering) with increasing AGDD, and the strength of this effect helps account for much of the trend seen in Fig. 3.

Table 3.

The Effect of Change in the Weather Covariates on the Timing of Flowering. We shifted each covariate by increasing its value each day by 30% of the mean across that day's range of variation in the historical record. Here we report the +30% shift, but the –30% shift in these covariates would produce the opposite effect on timing of flowering ( Fig. S2 (473_madr-68-04-18_s01.docx)). Delays and advances are reported. Unmodeled effects are indicated as empty cells (“-”). See  Figs. S1 (473_madr-68-04-18_s01.docx) and  S2 (473_madr-68-04-18_s01.docx) for more detail on effect size.

Fig. 5.

Contrasting effects of moisture on flowering times across different species. Panels A–B illustrate the expected influence of rainless days on the distribution of flowering times in the two years representing the extreme ends of rainless days for each species, with all other covariates held at their daily means. Panels C–D show the expected influence of vapor pressure deficit (VPD) on the distribution of flowering times in the two years representing the extreme ends of VPD for each species, with all other covariates held at their daily means. These results are based on the medians of all coefficient values reported in  Fig. S1 (473_madr-68-04-18_s01.docx). We show relative time to flowering (as opposed to date) along the x-axis here so that effect size is more apparent. See Table 1 for start dates.

Discussion

Our study sought to advance our understanding of phenological changes of keystone Sonoran Desert species and their response to climate variability nearly thirty years after Bowers and Dimmitt (1994) laid the foundation for dryland plant phenology. Using a new Bayesian model that sharpens our analysis of a unique long-term phenological data set, we demonstrate that the flowering phenology of keystone species in the Sonoran Desert are generally advancing and tied to rapidly changing temperatures – specifically more rapidly accumulating growing degrees early in the year. These results are consistent with others who have documented an advance in onset of flowering in woody species within the Sonoran Desert (Bowers 2007; Crimmins et al. 2010).

A data set collected by a single observer over a period of seventeen years immediately prior to the start of our phenological observations provides another opportunity to corroborate our results. McGinnies (1983) recorded the onset, peak and end of flowering for five of our six focal species (all but A. deltoidea) from 1966 through 1982 in a location 15 miles east of our study site. Comparing expected mean times for onset of flowering over our study period (1982–2018) with those reported by McGinnies for the seventeen years prior (1966–1982), we find that since the early 1980's, the onset of flowering has advanced for every species by days to weeks (2–16 days in our dataset, with an average across all six species of ∼9 days).

Flowering phenology of each species was also affected by short-term (2-week) windows of rainfall and VPD. However, unlike each species' consistent responses to variation in higher AGDD (advance), responses to moisture-based variables diverged across species. For example, whereas an increase in VPD advances flowering in O. tesota, it delays flowering in P. microphylla ( Fig. S2 (473_madr-68-04-18_s01.docx)). Some of the divergent responses across species to moisture-based variables may be explained by variation in functional traits that render them more or less sensitive to changes in temperature and moisture. Despite differences in responses to moisture variables, years in which flowering was exceptionally early or late were consistent across these tree species. For example, for both O. tesota and P. microphylla, the latest onset of flowering was observed in 1983 (Fig. 3). This consistency is likely a result of a similar response to AGDD, which outweighs divergent responses to moisture-related variables.

In contrast to our findings of an approximately two-day advance per decade from 1982 through 2018 for P. microphylla, Crimmins et al. (2010) noted a delay in flowering from 1984 through 2003. These differences may reflect a difference in two populations of P. microphylla, with our study population in the drier Tucson Mountains, and the other population in the more mesic foothills of the nearby Santa Catalina Mountains. Advancing flowering in our drier site may help ensure water availability during cooler months when potential evapotranspiration is low. An alternative explanation is that the trend has changed direction with the rapid warming in the additional fifteen years since 2003.

In comparison to the tree species, E. farinosa and A. deltoidea, two small shrubs, showed greater sensitivity to other climate drivers, including PDO and rainless days, and also antecedent monsoonal rainfall for A. deltoidea. While these shrubs can effectively modulate their water-use efficiency (Driscoll et al. 2020), they may be more sensitive to changes in water availability than our other species. Although AGDD pushes E. farinosa towards earlier flowering, any trend in the onset of flowering may be countered by its ability to respond to short-term rainfall (Fig. 5A). In contrast to the other five species, the probability of flowering in E. farinosa does not exhibit a distinct peak, but rather gradually rises and hovers around a near constant probability for at least three months before gradually declining (Fig. 4).

Mode of pollination may also constrain flowering phenology. The probability of flowering in E. farinosa and A. deltoidea exhibits the opposite response to rainless days. In A. deltoidea the probability of flowering increases with increasing number of rainless days, dropping precipitously after significant rainfall events (Fig. 5B), perhaps an adaptation for this wind-pollinated species. The opposite pattern is seen in E. farinosa, an insect-pollinated species. Its probability of flowering increases after each rainfall event (Fig. 5A). More generally, many of the dominant species in our study are sparsely distributed across landscapes in the Sonoran Desert and rely on pollinators, especially native bees (Simpson and Neff 1987). Future climate change could disrupt plant-pollinator interactions thereby affecting the reproductive output of these plant species.

The Sonoran Desert, like many global drylands, has warming rates much higher than other regions due to low vegetation cover and low soil moisture that enhance increases in temperature and aridity (Huang et al. 2016). Climate projections show a large amount of warming and drying expected for this site, with an increase of around 6 °C in mean annual temperature and more than 1 kPa in mean annual VPD, compared to the 1980s. Although the Sonoran Desert has experienced a reduction in mean annual precipitation since the 1980s (Munson et al. 2012), future projections show that precipitation may increase by 280 mm or decrease by as much as 120 mm. Regardless of the directional change in precipitation, increases in temperature and VPD alone will drive down soil moisture, thereby creating additional stress in plant water availability (Grossiord et al. 2020). Flowering phenology may continue to respond to climate trends moving forward, with higher divergence up to a limit, and pollinator activity and other non-climatic factors may interact to constrain future phenology.

With increasing temperature, we can generally expect an increase in accumulated growing degree days and advancements in flowering for these keystone desert species in the future. However, our models suggest that advancements in flowering may be accentuated or countered by concurrent changes in VPD, freezing events, and rainfall, depending on the species. Importantly, our uncertainty in the effect of any given covariate value on mean flowering time increases as covariate values become more extreme, which may be the case under future climates. When covariate values are more extreme (i.e., further from the overall mean covariate value in either direction) a relatively small amount of uncertainty in model parameters propagates into much larger uncertainty in outcomes. In other words, as covariate values become more extreme, there is less confidence in model predictions, and accordingly, variance is higher. This artifact (higher uncertainty and greater spread in possible mean flowering times) can be seen in our predictions for certain years that had extreme covariate values (Fig. 3).

Multi-decadal observational records and the increasing availability of powerful modeling techniques have improved our knowledge of how moisture availability affects flowering phenology, but our understanding remains incomplete. A focus on VPD in flowering phenology is particularly new, but initial work shows that developing and maintaining flowers is water-costly, especially in desert regions, for two reasons. First, water use during flowering can be high if flowers transpire more water than leaves, including corollas with greater surface area than leaves (Lambrecht 2013; Berg et al. 2019). Second, water content in corollas can be higher in some flowers that have a greater pectin content to maintain flowers, a common pattern in arid regions (Teixido et al. 2019). Our results provide evidence of both advancing or delaying floral development in response to increasing VPD, which is a way to avoid flowering during a time with high moisture deficit. The relationship between atmospheric water demand and the timing of these events is likely complex, including a plant's change in water status and functional need to allocate water rather than a direct relationship with atmospheric demand (e.g., Reich and Borchert 1984). Regardless, including short-term rolling windows of both VPD and rainfall, in addition to temperature, improved our phenological models, suggesting moisture availability in all its various forms, from dry soils to a thirsty atmosphere, is relevant to phenological processes and demands greater mechanistic understanding. This issue is especially true now that atmospheric demand is continuing to increase under climate change.

Higher temperature and growing degree days are associated with earlier pollinator emergence and metabolic and physiological processes that induce flowering (Hegland et al. 2009). Although the physiological processes that govern development in both insects and plants are regulated by temperature, this does not ensure phenological shifts in response to increasing temperatures will be similar in size for insects and their host plants (Forrest and Thomson 2011). From a pollinator's perspective (e.g., a Monarch Butterfly), this can result in reduced food availability during reproduction or migrations. While we did not monitor pollinator activity in our study, subsequent studies can help reveal to what degree shared phenological trends are driven by attracting a similar suite of pollinators and whether or not advances in flowering will impact the probability of successful pollination. Species likely to be most affected by changes in pollinator availability are obligate outcrossers, including Encelia farinosa (Ehleringer and Clark 1988), though plant population reductions may not be immediately realized for other long-lived Sonoran Desert species.

We document a trend of advanced onset of flowering in five of our six focal species and divergent responses to moisture-based variables. Much work remains to understand the implications of this trend and other phenological changes (e.g., duration of flowering) for the fitness and demography of our focal species and the implications for the broader ecological community. Species that inhabit highly variable environments, such as the Sonoran Desert might be expected to exhibit greater adaptive plasticity. Within just these six species, we observed a remarkable diversity of strategies to cope with limited and highly variable moisture availability. Such species may be more likely to survive in novel environmental conditions created by climate change.

Conclusions

Bowers and Dimmitt (1994) were among the first to quantitatively demonstrate changes in the phenology of Sonoran Desert woody species. Nearly three decades later, botanists at the Arizona-Sonora Desert Museum continue to collect phenological data with the same methodology on the same plant populations resulting in one of the longest records of arid plant species phenology in the world. We applied a new Bayesian model to an additional two decades of monitoring to take a more nuanced look at the phenology of these plant species. The new modeling approach we used is an improvement over simple correlative approaches and can advance future phenological studies by examining effects on daily probabilities of flowering, and by extension, the full distribution of flowering onset times. Importantly, improvements in modeling allow us to consider climate drivers at multiple temporal scales, from daily to antecedent. Future work could improve on this study by modeling intensity of response (the proportion of plants flowering) rather than simple presence / absence alone.

We provide evidence of a trend towards earlier flowering (1–4 days per decade) for five of the six species analyzed, likely linked to rising temperatures. Although dominant Sonoran Desert species similarly advanced in our study in response to greater AGDD, we found – like Bowers and Dimmitt (1994) – that measures of moisture availability were also important determinants of flowering phenology for some species that either countered or enhanced the advancing trends. The dominant species we evaluated showed either no effect or contrasting phenological responses to the same metrics of moisture availability, suggesting individualistic changes in phenology and the potential for divergent plant community flowering patterns under future climate change.