Herbarium specimens provide a critical source of phenological data that can be used to identify the direct and indirect drivers of variation in flowering date within and among species. Specimen-based phenological research in California has been accelerated by digitization efforts such as the California Phenology Network, which has scored and archived the phenological status of over 1.4 million specimens to date. Using this new data source in the Consortium of California Herbaria's CCH2 data portal, we obtained data from 993 specimens of the iconic California Poppy, Eschscholzia californica Cham., along with climate data representing all collection sites. Our goal was to determine how long-term and interannual climate variation affect flowering dates, and whether the magnitude of phenological sensitivity to climate varies across the species' range. We found that specimens collected from chronically warm or dry sites flowered relatively early, and flowering date was more sensitive to long-term mean temperature than to long-term mean precipitation. Independent of these effects of long-term conditions, flowering date in E. californica was sensitive to interannual variation in seasonal precipitation, but the direction of this effect depended on the season in which the precipitation occurred. Specimens sampled from sites experiencing warmer-than-average springs in the year of collection flowered 2.7–3.3 days earlier for every 1°C increase in spring temperature relative to long-term mean spring temperature. The magnitude of these effects, however, varied across the range of E. californica, with greater sensitivity to temperature in relatively cooler regions and no discernible sensitivity in relatively warm regions. Consistently, California Poppies exhibited significant phenological advancement over the last 120 years, but this advancement was restricted to the cooler portions of its range. Our results provide one of the first accounts of intraspecific variation in both phenological sensitivity to climate and the magnitude of phenological shifts over time, and demonstrate that, for a single species, location- or population-specific estimates of phenological sensitivity or of temporal trends in phenology might not accurately predict phenological responses to climate change in other locations throughout its range. In this study, we highlight the utility and promise of herbarium specimens for addressing novel questions about the phenological responses of plants to climate trends.
Shifts in plant phenological events (e.g., flowering or fruiting times) are well-documented effects of climate change across the globe (Parmesan and Yohe 2003; Fitchett et al. 2015; Scheffers et al. 2016; Piao et al. 2019). Because many organisms and ecological processes depend on the timing of plant phenological events, shifts in these events may affect species interactions and coexistence (Hegland et al. 2009; Kharouba et al. 2018; Rudolf 2019), ecosystem functioning (Richardson et al. 2010), and plant fitness (Mohan et al. 2019), especially since the phenological sensitivities of individual plants or populations may differ from those of the species with which they interact, even indirectly (Thackeray et al. 2016). For example, when elderberries and sockeye salmon became phenologically synchronized in Alaska due to differing responses to warmer spring temperatures, Kodiak brown bears consumed more elderberries and fewer salmon, altering the mortality rates of salmon and consumption of elderberry fruits (Deacy et al. 2017).
Differences among plant species in phenological sensitivity to climate have been widely reported (Wolkovich et al. 2012; Delgado et al. 2020; Parker 2021). However, the influence of climate on phenological events may also vary within species (Song et al. 2020), further complicating the prediction of the phenological effects of upcoming climate change, as these effects may be highly population- or region-specific. For several reasons, different parts of a widespread species' range might exhibit different degrees of phenological sensitivity to local, long-term climatic conditions or to inter-annual variation in climate. For example, with respect to temperature-sensitivity, in portions of a species' range where growth and reproduction are limited by cooler temperatures, we may expect populations to evolve strong sensitivity to inter-annual variation in temperature, resulting in earlier flowering in warmer-than-average years relative to “normal” years. Such phenotypic plasticity may be favored by natural selection because it enables individual plants to take advantage of the longer growing period that characterizes relatively warm years. By contrast, in portions of a species' range characterized by chronically warm, frost-free conditions, populations may evolve to flower in response to cues other than temperature, simply because local temperature might not be correlated with any particular risk through which natural selection would operate. Alternatively, if the relationship between local temperature or rainfall and flowering date is non-linear within a given species, then the phenological sensitivity to variation in these parameters would also necessarily change across temperature and/or precipitation gradients. In this case, even in the absence of genetically-based variation in phenological sensitivity to climate across the species’ range, we would observe regional variation in quantitative estimates of such sensitivity. In addition, if biotic factors that may influence the evolution of flowering time, such as the timing of pollinator availability or the timing of floral predation, covary with climate in some portions of the range but not others, then phenological sensitivity to climate would also appear to vary geographically. Whatever process drives spatial variation in phenological sensitivity to climate, it is important to document it, as accurate, region-specific forecasts of climate-driven phenological change depend upon both accurate predictions of the magnitude and direction of upcoming climate change and accurate estimates of phenological sensitivity to changes in temperature and precipitation.
To date, intraspecific variation in phenological sensitivity has been assessed in only a few taxa (Wang et al. 2015; Matthews and Mazer 2016; Prevéy et al. 2017; Park et al. 2018; Song et al. 2020), and the processes generating such variation are poorly understood. In the current study, we leverage the broad spatiotemporal coverage of digitized herbarium specimens of California Poppy (Eschscholzia californica Cham.) to detect and to interpret geographic variation in phenological sensitivity to three sources of climatic variation: geographic variation in long-term mean climatic conditions, inter-annual variation in temperature and rainfall, and seasonal variation in temperature and rainfall.
Our understanding of phenological shifts in plants has been derived largely from recurrent, location-specific observational records of living plants (Wolkovich et al. 2012, and references therein). These studies have advanced our understanding of the phenological responses of many species in a relatively small number of locations, but extrapolating to other regions or taxa is complicated by the magnitude of variation observed among taxa (Calinger et al. 2013; Park et al. 2018; but see Mazer et al. 2021). Furthermore, because observational phenological studies have predominantly estimated species-specific phenological sensitivities to climate using records from few locations, it is typically unclear whether they accurately represent the phenological sensitivity of a species in other locations within its range. To further understand the factors influencing phenological shifts within and among taxa, large quantities of data are required for many taxa representing different functional groups, with extensive coverage across time and space. Only in recent years have such data become available through herbarium digitization efforts (e.g., Yost et al. 2018) and observational data-collecting networks (e.g., the USA National Phenology Network, usanpn.org), enabling investigators to evaluate the causes and consequences of phenological shifts in rare or understudied plant communities (e.g., Munson and Sher 2015), across a more diverse set of taxa (e.g., Park and Mazer 2018), and across populations within individual species (e.g., Song et al. 2020).
Herbarium specimens, in particular, have emerged as a key source of phenological data, empowering studies across large spatiotemporal scales and for diverse and understudied plant, lichen, algal, and fungal taxa (Davis et al. 2015; Willis et al. 2017; Jones and Daehler 2018). Herbarium specimens have been collected across the world for several hundred years as vouchers for studies in taxonomy, biogeography, ecology, and other fields; their immense value for phenological research is a serendipitous byproduct of this foundational work. Herbarium specimens provide a “snapshot” of the phenological status of a given plant in a specific location on a certain date. In aggregate, analyses of specimen-based data have revealed how plants in different regions and environments (Park et al. 2018, Song et al. 2020), taxonomic groups (Calinger et al. 2013; Park et al. 2018), and functional groups (Calinger et al. 2013; Bertin 2015) respond phenologically to climate. In some cases, this work has identified the climate variables and seasons that are most influential for the studied species (Matthews and Mazer 2016; Park and Mazer 2018). Phenological estimates derived from herbarium specimens have demonstrated high reliability and congruence with those from other sources (Miller-Rushing et al. 2006; Panchen et al. 2012; Davis et al. 2015) despite known sampling biases (Daru et al. 2017).
Specimen-based phenological research, and specimen-based research in nearly all fields, has been greatly accelerated by digitization: imaging, transcribing label data, and georeferencing the collection locations of herbarium specimens (Willis et al. 2017). The National Science Foundation (NSF) has funded the mass digitization of natural history collections in the United States through the Advancing Digitization of Biodiversity Collections (ADBC) program, established in 2010. Over the last ten years, more than 125 million specimen records have been made freely available electronically (e.g., via iDigBio, https://www.idigbio.org/). These activities not only allow greater access to basic specimen data, but they also allow researchers to harvest phenological data on a larger scale and explore new avenues of research (Soltis 2017) such as those involving machine learning (Pearson et al. 2020).
Herbarium specimen digitization in California has a relatively long history, owing in part to the Consortium of California Herbaria (CCH), founded in 2003. Digitization activities have been recently accelerated by the California Phenology Thematic Collections Network (CAP TCN; https://www.capturingcaliforniasflowers.org/), funded by the NSF in 2018. This four-year project funds the digitization of nearly one million California vascular plant specimens (Yost et al. 2019). As a result, the project's new data portal, CCH2 ( https://cch2.org/) currently contains phenological data for over 1.4 million publicly accessible specimen records. By making these data available, the CAP TCN aims to facilitate phenological research that will improve our understanding of future changes to the California flora and help determine how potential threats to California biodiversity, ecology, and agriculture may be mitigated.
In this study, we used specimen occurrence data and newly created phenological data for 993 specimens of the California Poppy (Eschscholzia californica), the state flower of California, available in CCH2 to examine the effects of climate on the timing of flowering across the range of this iconic species. Using these data, we (1) identified the seasonal windows during which temperature and precipitation best predict flowering date of E. californica, (2) determined the magnitude and direction of the effects on flowering date of long-term temperature and precipitation, (3) assessed the effects on flowering date of interannual variation in seasonal temperature and precipitation, (4) detected whether temperature and precipitation have changed over time throughout the range of E. californica, and (5) examined whether the phenological sensitivity of this species to interannual climate variation or rates of phenological change over time vary among regions within its range. Additionally, although most studies have reported stronger phenological effects of interannual variation in temperature than of precipitation (Abu-Asab et al. 2001; Gordo and Sanz 2010; Hart et al. 2014; but see Matthews and Mazer 2016; Munson and Long 2017; Wolf et al. 2017), it is unclear whether the relative influence of these climatic factors is the same for species in water-limited regions, such as the native range of E. californica, or whether their relative importance is stable across a species' range. Therefore, we also compared the relative effects of interannual variation in temperature and interannual variation in precipitation on flowering date within the range of E. californica. These analyses provide one of the first assessments of intraspecific variation in phenological sensitivity to climate in Mediterranean environments and underscore the utility and promise of CCH2's herbarium specimen data for addressing novel questions concerning plant phenological responses to climate change.
Methods
Study System
The California Poppy, Eschscholzia californica, is a common and widespread annual or perennial herb that is native to the west coast of North America and has been introduced as a cultivated wildflower across the continent. California Poppies inhabit a wide range of habitats and climates but are most common in open grasslands, hillsides, coastal bluffs, and road cuts, where they co-occur with grasses and other wildflowers. The species is highly morphologically variable and has been variously divided into several taxa, yet all share the reproductive structures characteristic of the taxon: smooth, pointed buds enclosed by a fused, cap-like calyx; flowers with four bright yellow or orange petals; and narrow, cylindrical fruits arising from a swollen, obconic receptacle (Hannan and Clark 2012). An herbarium specimen displaying these features is shown in Fig. 1. We selected E. californica for this study because of its distinct and conspicuous reproductive features (which make the phenological status of specimens relatively easy to score), its wide distribution across California, and its iconic status in the state. Furthermore, emerging research suggests that phenological traits of E. californica are correlated with fitness and vary among populations occupying different climates (Ryan and Cleland 2019).
Dataset Preparation
We visualized all available and georeferenced California specimen records for E. californica using the Map Search feature of the Consortium of California Herbaria CCH2 data portal (2020; https://cch2.org) and identified regions of the state with relatively low sampling. We then georeferenced additional specimen records from poorly represented counties using GEOLocate ( https://www.geo-locate.org/) in CCH2 according to CAP TCN protocols (Wieczorek et al. 2012; Pearson 2020; Zermoglio et al. 2020). As a result, 2199 of the 4197 U.S. E. californica specimens in CCH2 were georeferenced and thus usable for our analyses. Taxonomic synonyms of E. californica (see Appendix S1 (343_madr-68-04-11_s01.docx)) were included in this search and treated as E. californica specimens for all analyses. We then used the Trait Coding from Images tool in CCH2 to score the phenological status of all E. californica specimens that had precise date and location information but that had not previously been scored by other efforts (see ‘Phenological data’ section below). Excluding specimens without precise date information reduced the dataset to 2101 specimens.
We downloaded the complete specimen records for the 2101 dated and georeferenced herbarium specimens of E. californica in CCH2 as of November 10, 2020. To download phenological data along with other available specimen data, we selected the “include Occurrence Trait Attributes” option when downloading the records from the CCH2 public search. This download results in a Darwin Core “MeasurementOrFact” extension file ( https://dwc.tdwg.org/terms/#measurementorfact) that is separate from the main occurrences table. Each line in the table corresponds to one phenological attribute (e.g., Open Flower present) that then corresponds to a specimen, such that several lines in the table may correspond to the same specimen. To merge the phenological data with the occurrences data in a more intuitive format, we used custom R code (R Core Team, R Foundation for Statistical Computing, Vienna, Austria) developed by the California Phenology Network ( http://doi.org/10.5281/zenodo.4298817).
We used the tidyverse package in R (Wickham et al. 2019) to filter the merged specimen and phenological data, excluding all specimens that lacked a phenological score (see ‘Phenological data’). We also removed specimens collected after day of year (DOY) 250 (September 7th) because we were primarily interested in factors influencing spring and summer flowering dates, and we suspected that flowering specimens collected after DOY 250 may have been either erroneously scored or collected from sites at which the plants were responding to disturbance, cultivation, or other non-seasonal cues. The resulting dataset consisted of 993 specimens.
Phenological data. The phenological scores of angiosperm specimens in CCH2 are coded according to the first- and second-order questions outlined by Yost et al. (2018). The possible first-order scores are “reproductive,” “sterile,” or “not scorable.” Second-order phenological scores reflect the presence or absence of four reproductive structures: unopen flowers, open flowers, senesced flowers, or fruits. Specimens with second-order scores that indicate the presence of any reproductive structure are automatically assigned a first-order score of “reproductive” as well.
Phenological scores in the Trait Attributes Table in CCH2 have primarily been produced by one of two methods: (1) visual inspection of specimens or specimen images or (2) text mining of label data. Visual inspection and scoring of specimen images method (1) is facilitated by the Trait Coding from Images tool in CCH2, which presents the user with images of specimens from any user-selected taxon. The user then selects the appropriate score(s) (i.e., presence or absence of reproductive structures) from a list of options displayed next to the specimen image. Phenological scores for E. californica specimens produced using this method were created by trained student technicians and staff.
For method (2), the text mining tool in CCH2 was used to view textual data from the labels of herbarium specimens and, when possible, to score the stated or implied presence of reproductive structures, as indicated by the text. For example, specimens with label text including “flowers red” would be scored as “open flower present.” Descriptions were interpreted conservatively; for instance, if a description included the phrase “red-flowered annual,” no phenological score was assigned because it was unclear whether the note refers to the species in general or the specimen itself. A screenshot of this tool is shown in Appendix S2 (343_madr-68-04-11_s02.docx).
We restricted our dataset to include only specimens scored as “open flowers present.” Although finer-scale methods of scoring the phenological status of herbarium specimens have been devised (Love et al. 2019; Pearson 2019a), prior research indicates that the coarse, presence/absence scoring schema is effective for identifying climatic trends given the large size of our dataset (Pearson 2019a). The resulting dataset consisted of 993 specimens, collected from 1901 to 2019, with high geographic variation in sampling location (Fig. 2).
Fig. 2.
Collection sites of the 993 E. californica specimens in our cleaned dataset. The shading gradient indicates long-term mean annual temperature (MAT) from 1901 to 2019 (see ‘Climate data’ for calculation of long-term MAT).
Climate data. We calculated long-term temperature and precipitation conditions in each collection site, and temperature and precipitation “anomalies” —the difference between climatic conditions in the year of specimen collection (as well as the previous year) and long-term mean temperature and precipitation for that site. To do so, we obtained monthly time series of mean temperatures and cumulative precipitation (hereafter ‘Tave’ and ‘PPT’, respectively) from January of 1901 to December of 2019 for each site of specimen collection in our sample using ClimateNA v.6.30 (Wang et al. 2016). ClimateNA is a software package that interpolates gridded data from PRISM (PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu) at a 4 km × 4 km resolution to produce scale-free climatic estimates adjusted for elevation. We calculated long-term climatic conditions by averaging mean annual temperatures (MAT) and cumulative annual precipitation (MAP) at each site of specimen collection throughout the study period (1901 to 2019). Similarly, for each collection site, we calculated long-term Tave and PPT for each calendar month (Jan–Dec) by averaging observed conditions for a given month between 1901 and 2019 (the full range of available climate data). We calculated longterm means over the full time range of the data to account for all climatic variability in our estimates of mean climate, but using a shorter, pre-warming window (1940–1970) generated nearly identical estimates. This yielded 120-year averages for MAT and MAP (hereafter ‘long-term MAT’ and ‘long-term MAP’, respectively), and 120-year averages for Tave and PPT for every calendar month (hereafter ‘longterm Tave or PPT’) in each collection site. For each specimen in our sample, we calculated monthly climate deviations (hereafter ‘anomalies’) from long-term conditions in its year and site of collection, and in the previous year, as the difference between observed Tave or PPT and long-term conditions for that month. This yielded 24 monthly climate anomalies (12 monthly anomalies in the year of collection, and 12 anomalies in the previous year) for both Tave and PPT for each specimen in our sample. No monthly anomaly was significantly correlated with its corresponding long-term mean (|r| < 0.31 for all months), indicating that multicollinearity did not affect our results.
Locations of specimen collection varied widely in long-term MAP conditions, with mean annual cumulative precipitation ranging from 220 mm to 1422 mm for the central 95% of the long-term MAP distribution. Consequently, the biological effects of net PPT anomalies of equal magnitude are likely to vary among sites of specimen collection because their magnitude in proportion to local mean precipitation conditions may be vastly different. To account for this, for each collection site, we calculated precipitation anomalies in each month in the year of collection, and the previous year, proportional to long-term MAP for that month. The resulting proportional PPT anomalies are expressed as a fraction of long-term MAP; a proportional PPT anomaly equal to 1 is of magnitude equal to the longterm MAP for that seasonal window in that site, whereas a value of 0.5 would correspond to an anomaly of 50% of the magnitude of long-term MAP for that period.
Analyses
Identifying phenologically influential seasonal windows. We identified the seasons during which PPT and Tave influenced flowering date most strongly using a sliding window analysis, separately estimating phenological sensitivity of flowering DOY (henceforth “flowering date”) to these climate variables in multiple time periods while controlling for variation in long-term MAP and in long-term MAT among sites. For each specimen in our sample, we used Tave and proportional PPT anomalies for each month (from September of the year prior to collection to May of its year of collection) as predictor variables. To obtain anomalies for broader time windows, we averaged monthly anomalies for successive 2-month (Sep–Oct, Oct–Nov, Nov–Dec, etc.) and 3-month periods (Sep–Oct–Nov, Oct–Nov–Dec, etc.). Then, we estimated the effects of PPT and Tave anomalies on flowering date in each time period by fitting a multiple regression for each time window and each climate variable. For each regression, we used the collection date of each specimen in year i and site j as a response (DOYij), long-term MAP and long-term MAT in each site j as controls (Long term MAPj and Long term MATj, respectively), and either Tave or proportional PPT anomalies in year j and seasonal window k (1-month, 2-month, or 3-month interval) as a predictor (Tave anomalyjk or PPT anomalyjk) (Equation 1).
Equation 1
Among time-windows for which climate anomalies produced statistically significant coefficients, those of greatest magnitude were interpreted as affecting flowering date most strongly in E. californica. Changes in the magnitude, direction, and significance of climatic parameters over consecutive time windows (e.g., shifts from negative to positive PPT coefficients from Fall to Spring months) were consistent among 1-month, 2-month, and 3-month windows. However, models including climate variables averaged over 2-month windows showed better fit (i.e., higher model R2 values) and lower standard errors in parameter estimation (data not shown). Consequently, we report results for time window analyses using 2-month periods.
Phenological sensitivity to long-term and interannual climate variation. To assess how variation in long-term climatic conditions and interannual climate variation influence the flowering date of E. californica, we evaluated the effects of Tave and proportional PPT anomalies and long-term MAP and long-term MAT conditions on flowering date among collection sites. We selected the 2-month windows for Tave and proportional PPT anomalies showing the greatest effects on flowering date. Some periods that showed significant effects for Tave and proportional PPT anomalies were moderately or highly correlated. Consequently, for both proportional PPT and Tave anomalies, we included more than one period only if (1) their correlation coefficients were less than or equal to 0.6, (2) inclusion of the period with an effect of lesser magnitude yielded a significant regression coefficient, and (3) inclusion of the additional period lowered the AIC value of the model by two or more.
To evaluate whether the effects of interannual variation in precipitation and temperature (i.e., anomalies) varied among parts of the range of E. californica with different long-term climatic conditions, we tested whether Tave and proportional PPT anomalies significantly interacted with long-term MAT and MAP across specimen collection sites. We fit multiple linear regression models including interaction terms for either long-term MAT or MAP and Tave or proportional PPT anomalies (i.e., MAT normals|Tave anomalies, MAT normals|PPT anomalies, MAP normals|Tave anomalies, and MAP normals|PPT anomalies for each 2-month window selected for inclusion in the final model). We only included one interaction term per regression model, and from these models we selected significant interaction terms for inclusion in the final model reported in this study. The final regression included flowering date for each specimen as a response (DOYij), long-term MAT and long-term MAP (long-term MATj and long-term MATj), April–May Tave anomalies (Apr–May Tave ij), February–March proportional PPT anomalies (Feb–Mar PPTij), October–November proportional PPT anomalies in the year (Oct–Nov PPTij), and the interaction between April–May Tave anomalies and long-term MAT (long-term MATj x Apr–May Tave ij) as predictors (Equation 2).
Equation 2
All predictors were mean-centered for easier biological interpretation of the intercept and main effects of the model. To compare the magnitude of predictor effects relative to their scale of variation, we fit the same model, scaling all predictors to a standard deviation of 1. By doing so, regression coefficients can be interpreted as the expected change in the response (flowering date as DOY) for a change of one standard deviation (SD) in each predictor. Variance inflation factors (VIFs) for all predictors were lower than 1.34, suggesting collinearity was unlikely to cause estimation problems.
Temporal trends in climate. We quantified temporal trends in the climate variables most strongly influencing flowering date in E. californica using generalized additive models implemented using the ‘mgcv 1.8-31’ package in R (GAMs; Wood, S. and M. S. Wood. version 1, 29). These models used Tave and proportional PPT anomalies as response variables (see Equation 2), year as a predictor, and longterm MAT and long-term MAP as covariates controlling for potential geographic variation in the locations of specimen collection over time. To assess whether the degree of climatic change varied throughout the range of E. californica, we repeated the analysis using linear regressions and included interaction terms between year and long-term MAT and long-term MAP. We found no significant interactions and report results only from the GAMs.
Temporal trends in phenology. We used multiple linear regression to assess whether flowering dates of E. californica have changed over time, and to evaluate whether rates of phenological change have varied among regions with different long-term climatic conditions within the species' range. Changes in the regions where specimens are sampled over time could lead to spurious temporal trends in phenology (e.g., recently collected specimens may come from warmer areas where flowering occurs earlier, on average, even if flowering dates have not changed over time); therefore, we included long-term MAP and long-term MAT as controls. To evaluate whether the magnitude of phenological change varies throughout the range of E. californica, we included interaction terms of year with long-term MAT and year with long-term MAP in each collection site, retaining only significant interaction terms in the final model. The resulting regression included flowering dates as a response (DOYij), year as a predictor (yeari), long-term MAT and MAP conditions as controls (Long-term MATj and Long-term MAPj), and an interaction term between long-term MAT and year of collection (Long-term MATj × yeari).
Equation 3
We tested for temporal and spatial autocorrelation in the residuals of all models using spatial and temporal distance matrix regressions, but found no evidence for either (R2 ≤ 0.001 and P ≥ 0.07 for all tests). Variance inflation factors (VIFs) for all predictors were lower than 1.29. All analyses were conducted using R v4.0.3 (R Core Team, R Foundation for Statistical Computing, Vienna, Austria). Interaction plots were visualized using the visreg package v2.7.0 (Breheny and Burchett 2017) and ggplot2 v3.3.2 (Wickham 2016).
Data Availability
The authors have made the cleaned specimen dataset and analysis code available via a public, online repository ( http://doi.org/10.5281/zenodo.4383254).
Results
Seasonal Phenological Sensitivities
The magnitude of phenological sensitivity of E. californica to temperature gradually increased from the fall of the year before flowering to the spring of the flowering year (Fig. 3a), with higher temperature anomalies associated with earlier flowering dates. We detected significant negative effects of temperature anomalies on flowering date (as DOY) for periods spanning the end of winter and most of the spring in the year of flowering, indicating earlier flowering in warmer-than-average years by 2.7–3.3 days/°C (Fig. 3a).
Fig. 3.
Sensitivity of flowering phenology to interannual variation in (a) mean temperature anomalies and (b) mean cumulative precipitation anomalies (proportional to mean precipitation in that period within each site) for different 2-month windows prior to the mean month of flowering of 993 herbarium specimens of E. californica. Asterisks represent significant phenological sensitivities to temperature or proportional precipitation anomalies for the focal 2-month period at α = 0.05 (*), α = 0.01 (**), and α = 0.001 (***).
We detected significant sensitivity to cumulative precipitation anomalies (proportional to long-term cumulative precipitation in each collection site) in both the fall of the year prior to flowering and the spring of the flowering year (Fig. 3b); however, proportional precipitation anomalies during these periods had opposite effects on flowering date, with greater-than-average September–October and October–November precipitation leading to advances in flowering date (–3.1 and –4.3 days per proportional PPT anomaly, respectively), and greater-than-average February-March and March-April precipitation resulting in delays in flowering date (6.8 days and 4.4 days per proportional PPT anomaly, respectively).
Phenological Sensitivity to Long-term and Interannual Climate Variation
Our model (Equation 2) explained 21% of the variation in flowering date among specimens (Table 1). We found significant relationships between longterm climatic conditions among sites of specimen collection and flowering date, with geographic gradients across increasing MAP associated with later flowering dates (0.8 days/100 mm; Fig. 4a), and geographic gradients across increasing MAT associated with earlier flowering dates (–7.2 days/°C; Fig. 4b). Among collection sites represented in our data, variation in standardized long-term MAT was associated with greater changes in flowering date than variation in standardized long-term MAP; changes in long-term MAT of one SD (2.0 °C) among collection sites were associated with average advances in flowering date of 15 days, while changes in long-term MAP of one SD (305 mL) were associated with average delays in flowering date of 2.6 days (Table 1, Figs. 4a and b).
Table 1.
Results of Linear Regression of Flowering Date and Climate Variables. Output of linear regression of day of year of flowering (DOY) vs. long-term mean cumulative precipitation (MAP) and mean annual temperature (MAT) at each site of collection, and cumulative precipitation and mean temperature anomalies for various 2-month periods in the year and site of specimen collection (See Equation 2 in Methods Section; adjusted R2 = 0.21, df = 986, F = 45.9, P < 0.001). All predictors in the model were mean-centered and standardized to a SD of one to facilitate comparison of effect sizes among predictor variables. Negative coefficients (Coef.) indicate advancement in DOY (earlier flowering), while positive coefficients indicate delays in DOY (later flowering).
Fig. 4.
Relationship between date of flowering (DOY) among 993 herbarium specimens of E. californica and (a) long-term mean MAP at the site of specimen collection (b) long-term mean MAT at the site of specimen collection (c) anomalies in cumulative February-March precipitation in the year of collection proportional to long-term averages for the site of collection, and (d) anomalies in cumulative October-November precipitation in the year of collection proportional to longterm averages for the site of collection.
The opposing effects on flowering date of proportional PPT anomalies for October–November and February–March revealed through sliding-window analyses remained consistent in magnitude and direction when both proportional PPT anomalies were included in the same model (see Equation 2). October–November and February–March proportional PPT anomalies had effects on flowering date of similar magnitude (but opposite in direction) per one SD of change in each variable (Table 1, Fig. 4c and d).
We detected a significant interaction between longterm MAT among sites of specimen collection and April–May Tave anomalies, indicating that the effects of interannual temperature variation on flowering date vary with long-term temperature conditions (Fig. 5). A difference of 2.4°C in long-term MAT relative to the mean among sites was associated with changes of ±3.0 days/°C in the effects of April–May temperature anomalies on flowering date (SE = ±0.96 days/°C). Accordingly, the phenological sensitivity to April–May temperature anomalies among E. californica specimens collected in relatively cold regions (2.4 °C lower than average, or 10th percentile) was –5.5 days/°C, whereas that of specimens collected in relatively warm regions (1.8 °C higher than average, or 90th percentile) was only –0.3 days/°C. The phenological sensitivity of E. californica to April–May temperatures and to October–November and February–March precipitation did not significantly interact with long-term MAP. Consequently, interannual temperature variation had greater effects on flowering date than interannual precipitation variation in colder-than-average sites, whereas precipitation had a greater effect on flowering date than temperature in warmer-than-average sites. For example, in relatively warm regions (1.8°C greater than average), an increase of 1.1°C (1 SD given the range of our data) in April–May temperatures was expected to result in a flowering advancement of only 0.4 days, whereas increases of 52% and 61% (1 SD given the range of our data) in fall and spring precipitation relative to long-term MAP in each collection site, which did not vary with either longterm temperature or precipitation, were associated with flowering advances of 3.7 days and delays of 4.0 days, respectively (Table 1). These results were consistent when using temperature anomalies calculated over other significant two month periods (February–March or March–April; results not shown).
Fig. 5.
Relationship between date of flowering (DOY) and mean April-May temperature anomalies across portions of the range of E. californica differing in long-term mean annual temperature (MAT). Red points and lines designate specimens in regions with warmer-than-average long-term MAT, green points and lines designate regions with approximately average long-term MAT, and blue points and lines designate specimens in regions with cooler-than-average long-term MAT.
Climate Change Over Time
We detected significant increases in Tave anomalies over time among the 993 collection locations represented in our sample, with an increase of approximately 1°C in mean April–May temperature from 1970 to 2019 (i.e., ∼ 0.2 °C/decade; P << 0.001; Fig. 6a). We detected a significant reduction in October–November proportional PPT anomalies over time among collection sites (P = 0.002; Fig. 6b). However, the reduction in proportional PPT anomalies between 1970 and 2019, the period best represented in our sample and that marks the onset of rapid climate change, was very small relative to the range of variation in the data; decreases in PPT anomalies proportional to long-term MAP in each site since 1970 were lower than 0.1, corresponding to a reduction of less than 10% of the mean October–November cumulative precipitation per site, whereas mean interannual variation was 0.52, or 52% of October–November long-term MAP in a site. We found no significant temporal trends in February–March proportional PPT anomalies among specimen collection sites (Fig. 6c).
Fig. 6.
Change over time in (a) April-May temperature anomaly (°C), (b) February-March cumulative precipitation anomaly (proportional to long-term averages for the site), and (c) October-November cumulative precipitation anomaly (proportional to long-term averages for the site) across 993 sites of specimen collection throughout the range of E. californica.
Temporal Trends in Flowering Date
Our model designed to detect temporal trends in E. californica flowering date (Equation 3) revealed an average advancement of flowering date of 0.8 days per decade (P = 0.048). However, we detected a significant interaction between long-term MAT and the effect of year, indicating that rates of phenological change over time across the range of E. californica depended on local long-term temperatures (Fig. 7, P = 0.002). A difference of 2.4°C in long-term MAT relative to the mean among sites was associated with changes of ±1.4 days per decade in the rate of flowering date change over time. Accordingly, specimens collected in colder-than-average regions (2.4 °C lower than average; 10th percentile) were predicted to advance their flowering phenology at a rate of 2.6 days per decade, whereas specimens collected in relatively warm regions (1.8 °C greater than average; 90th percentile) were predicted to delay flowering at a rate of only 0.4 days per decade.
Fig. 7.
Temporal trend in the flowering phenology of E. californica across portions of its range differing in long-term mean annual temperature conditions. Red points and lines designate specimens in regions with warmer-than-average long-term MAT, green points and lines designate regions with approximately average long-term MAT, and blue points and lines designate specimens in regions with cooler-than-average long-term MAT.
Discussion
Differences in long-term temperature and precipitation were associated with differences in flowering date, with later flowering at sites with relatively high precipitation, and earlier flowering at warmer-than-average sites. Furthermore, flowering date variation depended more strongly on long-term mean temperature than on long-term mean precipitation. Interannual variation (i.e., anomalies) in both temperature and precipitation influenced flowering date, but this sensitivity depended on long-term temperature at the collection site, with greater sensitivity to temperature anomalies in colder regions. Consistent with observed temporal warming trends across the range of E. californica, we detected an advancement of flowering date over time that mirrors the long-term temperature-dependence of phenological sensitivity to changes in temperature; the rate of change of flowering date in cooler-than-average sites has outpaced that of warmer-than-average sites. This intraspecific variation in phenological sensitivities across the range of E. californica, along with regional differences in projected climate change, may result in complex phenological changes for E. californica that may impact fitness and species interactions within herbaceous forb and grassland communities. The spatiotemporal coverage of our data, enabled by the digitization of herbarium specimens, provides a nuanced view of range-wide phenological trends that may inform future studies designed to distinguish the mechanisms and effects of phenological shifts.
Phenological Sensitivities of Flowering Date to Seasonal, Interannual Variation
Of the interannual seasonal temperature variables investigated in this study, late-winter and spring (Feb–May) temperature anomalies in the year of flowering had the strongest effects on flowering date of E. californica. Increases in temperature advanced flowering at a rate of 2.7–3.3 days per 1°C increase in spring temperature anomaly (Fig. 3a). Fall and winter temperature anomalies in the year prior to flowering did not significantly affect flowering date. This may reflect a short or nonexistent vernalization requirement in E. californica, which could be expected for plants in Mediterranean climates that are unlikely to experience very cold winters. The advancing effect of spring temperature anomalies on flowering date is consistent with estimates of flowering advancements in other climates and for other taxa (Miller-Rushing and Primack 2008; Willis et al. 2017; Pearson 2019b). Considering the warming trends that we observed across the range of E. californica (Fig. 6) and projections of future climate change (Bedsworth et al. 2018), flowering dates of E. californica are projected to continue to advance.
However, our analyses also revealed that phenological sensitivities to interannual temperature variation depended on long-term temperature at the site of collection; E. californica in cooler-than-average sites were more responsive to spring temperature than those in warmer-than-average sites. This site-dependence resulted in phenological shifts of several days per °C in colder-than-average locations, and in an apparent lack of phenological responsiveness to temperature in warmer-than-average locations. Few herbarium-based studies have quantified intraspecific variation in phenological sensitivity to climate and of phenological change over time across the range of a single species (but see Matthews and Mazer 2016; Munson and Long 2017; Song et al. 2020, 2021). Our results demonstrate that, for a single species, estimates of phenological sensitivity derived from data in a single location might not accurately predict phenological sensitivities to climate in other locations throughout the range of the species. Moreover, patterns of intraspecific variation may differ among taxa, making generalization difficult. For example, Song et al. (2020) discovered greater phenological sensitivity of a perennial herb in dry, warm climates, in contrast to our findings. The mechanisms underlying the opposite trend in E. californica are unknown but may be related to developmental limits; California Poppies in warmer-than-average locations may already be flowering as early as developmentally possible, and flowering any earlier might mean doing so at such a small size that the amount of stored resources is insufficient to ensure fruit maturation. Alternatively, flowering earlier could result in phenological mismatches with their pollinators. If either scenario applies, then poppies in warmer climates may face particular challenges with future climate change. If poppies are unable to shift flowering into cooler months or otherwise adapt to environmental stressors, extreme heat or drought during flowering and fruit development could increase mortality and decrease individual fitness by decreasing reproductive output (Cox. et al. 2021). In these circumstances, the distribution of E. californica may shift toward cooler regions at higher latitudes and elevations, similar to what has been observed and predicted for other plant taxa in California (Kelly and Goulden 2008; Ackerly et al. 2020).
Predicting flowering date of E. californica in response to climate conditions is further complicated by the opposing effects of proportional cumulative precipitation anomalies in different seasonal windows. Flowering date was significantly delayed in response to increased spring (Feb–Apr) proportional PPT anomalies, but advanced with increased fall (Oct–Nov) proportional PPT anomalies in the year prior to flowering (Fig. 3b). Previous research has discovered contrasting effects of seasonal temperatures on phenological events (Cook et al. 2012; Hart et al. 2014; Pearson 2019b); however, to our knowledge, differences in the direction of the effects of precipitation for different seasons have seldom been reported (but see Munson and Sher 2015). Precipitation has commonly been found to have only a delaying or non-significant effect on flowering date (Abu-Asab et al. 2001; Matthews and Mazer 2016; Pearson 2019b). In contrast, interannual variation in both fall and spring precipitation influence flowering date in E. californica. Increases in fall precipitation might spur earlier germination (for annual individuals) and accelerate vegetative growth of E. californica (for both annual and biennial individuals), ultimately advancing spring flowering dates. However, higher-than-average precipitation in the spring may decrease the activity of pollinators (Lawson and Rands 2019) and extend the duration of the growing season; consequently, E. californica might have evolved to delay reproductive onset in springs with relatively high precipitation, potentially resulting in greater vegetative growth and resource accumulation prior to flowering.
Relative Influence of Temperature Versus Precipitation on Flowering Date
Although our analyses found that interannual variation in precipitation affected flowering dates, long-term mean temperature was the strongest predictor of flowering date, strongly outweighing the effects of long-term mean precipitation (Fig. 4, Table 1). This finding is consistent with previous research in Mediterranean (Gordo and Sanz 2010) and other climates (Abu-Asab et al. 2001; Hart et al. 2014; Munson and Long 2017), yet differs from other studies in the western U.S. (Matthews and Mazer 2016). Given the extreme water limitation in some areas of the range of E. californica and the potential effects of this limitation on fitness and natural selection (e.g., Franks et al. 2007), an overall lower impact of precipitation on flowering dates is somewhat unexpected. It is possible that, rather than having strong effects on flowering date, precipitation affects the duration of flowering and, consequently, total flower, fruit, and seed production. The causal relationships between flowering date and long-term temperature and precipitation are difficult to identify in observational studies. The relationships reported here may reflect adaptation to local climate conditions or they may reflect evolutionary or plastic responses to other environmental factors that may co-vary with climate, such as nutrient availability, pollinator abundances, and exposure to herbivores. Our analyses cannot distinguish between adaptive and plastic phenological responses to climate variation.
Although long-term temperature was consistently more predictive of flowering date than long-term precipitation, the effects of temperature and proportional precipitation anomalies on flowering date depended on long-term temperature at the collection site. Temperature anomalies had a greater effect on flowering date in colder-than-average sites, whereas precipitation had a greater effect on flowering date in warmer-than-average sites. This indicates that water limitation may have played an important role in shaping the phenological sensitivities of E. californica to interannual climate variation, and water availability will likely shape the phenological impacts of future climate change. Given the species- and region-specific effects of climate on phenological events (Petrauski et al. 2019; Song et al. 2020), more research is needed to investigate the relative importance of temperature and precipitation on California plants.
Temporal Trends: the Combined Effects of Phenological Responses to Long-term and Interannual Variation in Climate
Together, the contrasting and context-specific phenological responses of E. californica to temperature and precipitation, coupled with significant climate change across the range of the species (Fig. 6), forecast complex future phenological changes for E. californica. Our analyses suggest that intraspecific variation in phenological sensitivity has already caused disparate temporal shifts in flowering dates across the range of E. californica; flowering dates in cooler-than-average sites have shifted an estimated 3.1 days earlier per decade between 1901 and 2019, while flowering dates in warmer-than-average sites have shifted only an estimated 0.1 days per decade over the same timeframe. It is important to note that, had we not estimated the independent and interacting effects of long-term means and anomalies of temperature and precipitation on flowering date, the magnitude of sensitivity to temperature anomalies and of shifts in flowering date would be, on average, much lower than the sensitivity we identified in colder-than-average regions of the range of E. californica. These results highlight the importance of quantifying intraspecific variation in phenological sensitivity throughout species' ranges, and suggest that estimates of phenological change over time based on data from single locations should be interpreted with caution. Likewise, our results show that the flowering date of E. californica is sensitive to interannual variation in both temperature and precipitation, and, for precipitation, variation in two different seasons (fall and spring) had opposing effects on flowering date. Because each of these variables had unique and sometimes opposing effects on flowering date, the net effect in some climatic conditions could be no change in flowering date even though the species is responsive to both temperature and precipitation. As a result, species that appear to be phenologically “insensitive” to past climatic conditions may not be so in all climatic conditions in the future (Hart et al. 2014), especially as changes in temperature and precipitation may have complex effects on soil moisture. This highlights the importance of large, spatiotemporally diverse phenological datasets such as those afforded by digitized herbarium specimens to tease apart intraspecific and spatiotemporal variation in the phenological effects of climate change.
Limitations and Strengths of Specimen-based Phenological Analyses
Specimen-based phenological data, and thus the analyses presented here, are subject to certain limitations. Herbarium specimens are often collected opportunistically rather than systematically and may therefore contain spatiotemporal and taxonomic biases, such as over-collection near populated places and roads (Dodd et al. 2016; Daru et al. 2017), which could cause under-representation of certain environmental conditions in our data. For example, 42% of specimens in our dataset were collected in the southernmost 10% of latitudes. Collectors may also preferentially collect individuals that are in flower even if the majority of the population is not flowering (Willis et al. 2017). These limitations are largely overcome by the size and provenance of the dataset, our dataset represents the combined collecting effort of over 600 collectors across 22 herbaria, and by data cleaning steps to minimize the effects of outliers that originate from collector bias. Specimen-based data, like many data sources, are correlational and can only hint at the mechanisms underlying spatiotemporal changes in phenological events. Further experimental studies are needed to distinguish how context-specific factors, such as functional traits and local microclimates, will affect phenological events in E. californica and other taxa across California's diverse landscape.
Our study focused on a single wildflower species in California, the California Poppy, yet this study represents a path forward to understanding the phenological effects of climate change on plant taxa in the biologically and geographically diverse California Floristic Province. With the availability of large phenological datasets from the digitization and subsequent phenological scoring of herbarium sheets, we are poised to investigate more deeply the species- and region-specific effects of climate on plant phenological events. Furthermore, these types of analyses will enable us to ask broader questions about change across the state, such as whether temperature is consistently a stronger predictor of phenological events in annual and perennial herbs across California, whether phenological sensitivity is phylogenetically conserved (Davies et al. 2013; Mazer et al. 2013), and whether phenological sensitivity differs among species occupying or adapted to unique geological conditions such as serpentine soils. The taxonomic, geographic, and temporal scope of digitized herbarium specimens has the potential to open a new frontier in the study of plant responses to climate change.
Acknowledgments
The authors gratefully acknowledge funding from the U.S. National Science Foundation (NSF) to JY (DBI-1802312) and to SJM (DBI-1802181), and from UCSB, which provided generous graduate support to TRP in the form of a Chancellor's Fellowship.
