To date, most herbarium-based studies of phenological sensitivity to climate and of climate-driven phenological shifts fall into two categories: detailed species-specific studies vs. multi-species investigations designed to explain inter-specific variation in sensitivity to climate and/or the magnitude and direction of their long-term phenological shifts. Few herbarium-based studies, however, have compared the phenological responses of closely related taxa to detect: (1) phenological divergence, which may result from selection for the avoidance of heterospecific pollen transfer or competition for pollinators, or (2) phenological similarity, which may result from phylogenetic niche conservatism, parallel or convergent adaptive evolution, or genetic constraints that prevent divergence. Here, we compare two widespread Clarkia species in California with respect to: the climates that they occupy; mean flowering date, controlling for local climate; the degree and direction of climate change to which they have been exposed over the last ~115 years; the sensitivity of flowering date to inter-annual and to long-term mean maximum spring temperature and annual precipitation across their ranges; and their phenological change over time. Specimens of C. cylindrica were sampled from sites that were chronically cooler and drier than those of C. unguiculata, although their climate envelopes broadly overlapped. Clarkia cylindrica flowers ~ 3.5 days earlier than C. unguiculata when controlling for the effects of local climatic conditions and for quantitative variation in the phenological status of specimens. However, the congeners did not differ in their sensitivities to the climatic variables examined here; cumulative annual precipitation delayed flowering and higher spring temperatures advanced flowering. In spite of significant spring warming over the sampling period, neither species exhibited a long-term phenological shift. Precipitation and spring temperature interacted to influence flowering date: the advancing effect on flowering date of high spring temperatures was greater in dry than in mesic regions, and the delaying effect of high precipitation was greater in warm than in cool regions. The similarities between these species in their phenological sensitivity and behavior are consistent with the interpretation that facilitation by pollinators and/or shared environmental conditions generate similar patterns of selection, or that limited genetic variation in flowering time prevents evolutionary divergence between these species.

This study comprised five steps. First, we sampled reproductive specimens of each species from throughout its range to record each specimen’s date and site of collection, as well as a quantitative estimate of its phenological status, which ranged from bearing only flower buds to bearing only ripe fruits. Second, we compared the climatic conditions occupied by each species to assess their habitat preferences and tested for species differences in mean flowering date independent of local, long-term climatic conditions. Third, we used linear models to detect the degree of climate change that each species experienced over the past 112 – 119 years across their sampled range (1900 – 2012 for C. cylindrica and 1892-2011 for C. unguiculata). Fourth, we constructed and tested phenoclimatic models to detect the effects of local mean maximum spring temperature (Spring Tmax) and annual precipitation (AP) on each species’ flowering date, which was estimated as the day of the year (DOY) on which a reproductive specimen was collected; these linear models also controlled statistically for each specimen’s phenological status (estimated using a quantitative index of its reproductive progression; Love et al. 2019).

The phenoclimatic models tested here provide measures of the sensitivity of flowering phenology to climatic conditions estimated at two temporal scales: decadal and year-of-collection.  Specifically, each site of specimen collection was characterized by its mean climatic conditions from 1921-2010 (i.e., long-term climate, or climate “normals”), and by the deviation between climatic conditions at the site in the year of collection and the site’s mean long-term climate conditions (i.e., climate anomalies due to inter-annual variation in temperature and annual precipitation [AP]). When both kinds of parameters are included as explanatory variables in linear models designed to predict the DOY of flowering, the sensitivity of DOY to climate normals reflects a combination of local adaptation and plastic responses to spatial variation in chronic climatic conditions, while phenological responsiveness to interannual variation is due primarily or wholly to plastic responses to short-term local conditions.  Other investigators have used this approach, including Mazer et al.’s (2020) study of seed size variation in Clarkia; Pearson et al.’s (2022) investigation of Eschscholzia californica (California Poppy; Papaveraceae), and Parker’s (2022) study of five species of Arctostaphylos (Ericaceae) and Ceanothus (Rhamnaceae).

In all models, we controlled for variation in the phenological status of specimens, which can confound the relationship between a specimen’s collection date and the actual date of flowering onset. These models were then used to compare the direction and magnitude of the two species’ phenological sensitivities to local temperature and precipitation normal and anomalies, and to test for interactions between climate variables that may have influenced flowering time. Finally, we tested for phenological shifts in estimated flowering date within each species during the ~115-year sampling period to determine whether long-term temporal trends in flowering phenology were consistent with each species’ sensitivity to inter-annual variation in temperature and precipitation, and with the degree of climate change that each species experienced.

Study species

Clarkia is well-studied genus of ~40 species of self-compatible, annual, herbaceous wildflowers native to the western U.S. (Lewis and Lewis, 1955). Wherever they occur, populations of Clarkia are among the last spring wildflowers to bloom (typically flowering in May or June), and they produce dense and showy floral displays. The two taxa selected for the current study – Clarkia cylindrica ssp. clavicarpa W. Davis (subgenus Peripetasma) and C. unguiculata Lindley (subgenus Phaeostoma) – inhabit open or disturbed habitats and roadsides in the foothills, grasslands, and oak/pine woodlands of the Coastal Ranges, Transverse Ranges, and Sierra Nevada in California. These taxa are adapted to a Mediterranean climate, although the sites sampled for the current study vary widely with respect to long-term conditions. Among sites, long-term MAP estimated from 1921-2010 ranges from 141 – 1377 mm, and mean spring Tmax for the same period ranges from 9.9-24.9 °C. Clarkia cylindrica has been described as “normally outcrossing” (Davis, 1970), and C. unguiculata is predominantly outcrossing, although populations differ in their outcrossing rates (Vasek, 1958; Hove et al., 2016; Ivey et al., 2016). Both species are restricted to California, are typically found at elevations below 1500m, and are diploid (n=9).

The two species have been found to co-occur regularly in the southern Sierra Nevada (Moeller, 2004), but they are also found with other congeners. Lewis and Lewis (1955) reported that, at the time of the publication of their 1955 monograph on Clarkia, C. unguiculata and C. cylindrica could be found in mixed colonies with, respectively, any of 17 or 10 other Clarkia species (including each other). So, the potential for these species to compete for pollinators or to facilitate the pollination of congeners is not restricted to their interactions with one another. In addition, Lewis and Lewis (1955) found these two species to be genetically incompatible; no hybrids were formed following interspecific pollination. This incompatibility means that, while divergent flowering times might help sympatric populations of C. unguiculata and C. cylindrica to avoid interspecific pollination and stigma clogging, phenological divergence is not needed to achieve reproductive isolation.

Both species are pollinated by a small number of specialist, solitary “Clarkia bees” (Lewis and Lewis, 1955; MacSwain et al., 1973) that are attracted to the conspicuously pigmented, flecked, and spotted cup-shaped flowers of C. cylindrica and the rotate, clawed-petalled flowers of C. unguiculata, both of which provide pollen and nectar rewards (Fig. 1). Clarkia cylindrica is commonly visited for its pollen by Andrena lewisorum Thorp, which has also been found to visit C. unguiculata (MacSwain et al., 1973). Other bee species that visit C. cylindrica and C. unguiculata (but not necessarily at the same locations, at similar frequencies per plant species, or with equivalent pollination efficacy) include: Andrena omninigra Viereck, Hesperapis regularis (Cresson), Megachile gravita Mitchell, M. pascoensis Mitchell, Diadasia angusticeps Timberlake, Melissodes clarkiae LaBerge, Synhalonia venusta carinata (Timberlake), Ceratina sequoiae Michener, Lasioglossum pullilabre (Vachal) and Apis mellifera Linnaeus (MacSwain et al., 1973).

Herbarium data

The two Clarkia species examined here are well represented in the holdings of the Consortium of California Herbaria. For the current study, we borrowed specimens of each species from the Jepson Herbarium (JEPS) and the University Herbarium (UC) at the University of California, Berkeley; Rancho Santa Ana Botanic Garden (RSA); the University of California, Riverside (UCR), and the Santa Barbara Botanic Garden (SBBG). Each specimen’s label was examined to extract its collection date (recorded as the day of year on which the specimen was collected [1-365 or 366 on leap years], year of collection, elevation (m), latitude and longitude. Specimens that lacked GPS coordinates were georeferenced using the label’s description of the collection site and the online utility, GEOLOCATE (http://www.museum.tulane.edu/geolocate/). Where elevation was not recorded on the label, we used the specimen’s latitude and longitude (with GEOLOCATE) to estimate its elevation. Only specimens that bore one or more individual plants that had begun to reproduce (i.e., bearing flower buds, open or spent flowers, expanding ovaries, or fruits) at the time of collection were included in this study. Individual plants that had not begun to reproduce were not scored.

Because herbarium specimens may have been collected at any time during an individual’s reproductive cycle, a specimens’ collection date (DOY) is not a precise proxy for date of onset of flowering. Moreover, DOY is generally positively correlated with an individual plant’s phenological status such that, under similar environmental conditions, individuals collected at early stages of reproduction (e.g., when bearing only closed flower buds) are represented by earlier DOYs than those collected at later stages (e.g., when bearing only ripe fruits). Because DOY is confounded with reproductive stage, predictive models that control for the phenological status of sampled plants when examining the relationship between DOY and climatic conditions explain a higher proportion of the variance in DOY than those that do not control for this source of variance (Love et al., 2019).

To provide a quantitative estimate of each specimen’s phenological status, we used a phenological index (PI) that ranges from 1 (for plants comprised entirely of buds) to 4 (for plants comprised of all fruits) (Love et al., 2019). For each individual plant specimen (including multiple plants when a given herbarium sheet contained more than one individual), we counted the numbers of buds (> 5mm in length), open flowers, wilted flowers or expanding ovaries, and fully developed fruits (full-sized and/or beginning to dehisce). Each organ type was assigned a value that reflected its developmental stage (buds=1; open flowers = 2; spent flowers or expanding ovaries = 3; fully developed fruits = 4) and used in the following equation:

where Px represents the proportion of reproductive organs in a given stage and i represents the value assigned to that class (e. g., buds have a value of 1). The PI is therefore the weighted average of the proportions of buds, open flowers, spent flowers or developing ovaries, and ripe fruits. For herbarium sheets on which multiple, complete individuals were pressed, we assigned the sheet the mean of the PI values of its component individuals. Partial individuals were not scored for their PI. For C. unguiculata, a total of 608 plants on 231 sheets were scored; for C. cylindrica, a total of 1042 plants on 306 sheets were scored.

Duplicate specimens were considered to be those that were: collected < 500 m away from the nearest specimen in both latitude and longitude; collected on the same day of the same year; and represented by the same mean annual temperature normal, mean spring maximum temperature normal, and annual precipitation normal (as extracted using the climate database, ClimateNA.) Two sites separated by latitudinal and longitudinal distances of 500 m would be a linear distance of ~707 m apart, and were usually associated with distinct climate variables due to a difference in slope, aspect, or elevation, which ClimateNA uses to estimate climatic parameters. Following the elimination of duplicate specimens, 226 sheets of C. unguiculata collected from 1902-2011 and 284 sheets of C. cylindrica collected from 1900-2012 were analyzed here.

Within each species, some specimens were retained even if they were collected < 707m from another collection site. For 85 specimens of Clarkia cylindrica, the distance to the nearest collection site was < 707m, but only 24 of these specimens (distributed in 10 groups of 2 or 3 specimens) were collected on the same day in the same year as a nearby specimen. Based on the climate data retrieved from ClimateNA, the 2-3 specimens in each of these groups were represented by different combinations of mean annual temperature, annual precipitation, and mean Spring Tmax (likely due to differences between them in slope, aspect, and/or elevation, for which ClimateNA interpolated distinct climate parameters), and so they were retained for analysis. Twenty-two of the 226 specimens of Clarkia unguiculata were collected < 707m from the next closest site. Of these specimens, two were collected on the same day in the same year as their nearest neighbor, but 40 meters apart in elevation.

The range of PI values recorded for the sheets of C. cylindrica was 1.0-3.77; the range for C. unguiculata was 1.0 – 3.95. In all analyses described below, the PI was log-transformed to improve normality.

Climate data

We evaluated climate variables commonly found to influence flowering date in other taxa (Anderson et al., 2012; Cleland et al., 2012; Park and Mazer, 2018; Berg et al., 2018): annual precipitation (AP; this includes cumulative rain and snow, with the latter converted to water-equivalents); mean maximum spring temperature (Spring Tmax, the mean maximum daily temperature from March-April), mean minimum spring temperature (Spring Tmin), and the number of frost-free days in winter and in spring (NFFD). For each collection site and climatic parameter, two values were extracted from ClimateNA (Wang et al., 2016), an application that assembles downscaled monthly climatic parameters (e.g., the mean of daily values) for a wide range of parameters recorded from 1901 onwards. First, we obtained the long-term mean value of each parameter from 1921-2010. Second, we extracted the climate conditions for the year of collection (YOC). For each collection site x year combination, we then calculated the deviation between the annual conditions in the YOC and long-term climate mean (hereafter referred to as the “normal”). The value of this deviation (referred to here as the “anomaly” for the focal variable) indicates whether, in the year of specimen collection, a given site was warmer- (or cooler) or drier- (or wetter) than its long-term average.

Exploratory linear models were constructed and tested to determine whether spring Tmax or spring Tmin better explained variation in DOY. Each of two models was tested separately on C. cylindrica and C. unguiculata, and then on the pooled data including both species. The first model included the following variables as predictors: Log(PI), AP normal, Spring Tmax normal, AP anomaly, and spring Tmax anomaly. The second included the same predictors but using the Spring Tmin normals and anomalies instead of Tmax. The models in which the predictor variables included Spring Tmax performed better (had higher R² values) than those that included Spring Tmin (Table A1). Because Spring Tmax and Spring Tmin are collinear among collection sites (Spring Tmax normal vs. Spring Tmin normal: r = 0.55, P < 0.0001, n=509; Spring Tmax anomaly vs. spring Tmin anomaly: r = 0.75, P < 0.0001, n=509), we did not include both in the same model.

In addition to examining the effects of Annual Precipitation, Spring Tmax, and Spring Tmin on flowering date, we used data from ClimateNA to calculate the anomalies for Winter Tmax, Winter Tmin, precipitation as snow (PAS), and the number of frost-free days (NFFD) in winter and spring in order to examine long-term temporal trends in climate over the ~115 years of specimen collection represented by the data analyzed here.

Analyses

Climatic conditions and geographic locations occupied by each species. To compare species with respect to the combinations of climatic and geographic conditions they occupy throughout their sampled ranges, we examined the bivariate space occupied by each species’ collection sites with respect to their Spring Tmax and AP normals, as well as their elevation, latitude and longitude. Linear models were then constructed to test for significant differences between species’ means with respect to each of the four focal climate variables (AP normal, Spring Tmax normal, AP anomalies, and Spring Tmax anomalies) while controlling for variation in the other three variables. In each of these models, the focal climate variable was the response variable and the other three climate variables and Species were included as fixed main effects. One-way ANOVAs were conducted to compare species’ means with respect to elevation.

Climate change during sampling period: Temporal change in climate anomalies. Analyzing each species separately over its ~115-year sampling period, we conducted simple regressions to test for temporal trends in each of our focal climate variables: AP, PAS, winter Tmin, spring Tmin, winter Tmax, spring Tmax, winter NFFD, and spring NFFD. In each regression, the anomaly for a given climate variable was included as the response variable and Year as the independent variable. In these analyses, positive slopes of the regression of the temperature- (or precipitation-) based anomalies on year mean that, as time progresses, the sampled sites are becoming warmer (or wetter) than average.

Factors influencing flowering date: species identity, spring Tmax and AP normals and anomalies, and phenological status. We constructed and tested a suite of linear models to detect significant differences between species in mean DOY and to measure the independent effects on DOY of the normals for AP and spring Tmax, of the AP and spring Tmax anomalies, and of the phenological status (estimated as log[PI]) of individuals. Because interactions between climate variables can influence the interpretation of their individual effects on DOY, we also sought evidence for significant interactions between each pair of the climate variables that were included as main effects. This analysis was conducted in several steps; in all models, DOY was the response variable.

The first model tested only for the following main fixed effects: Species, Log(PI), AP and spring Tmax normals, and AP and spring Tmax anomalies. We then tested each of 15 two-way interactions by adding one of the following interactions to the first model to test its contribution to variance in DOY: AP normal|Tmax normal, AP normal|AP anomaly, MAP normal|Tmax anomaly, Tmax normal|AP anomaly, Tmax normal|Tmax anomaly, AP anomaly|Tmax anomaly, LogPI|Species, AP normal|Species, spring Tmax normal|Species, AP anomaly|Species, Spring Tmax anomaly|Species, AP normal|Log(PI), Spring Tmax normal|Log(PI), AP anomaly|Log(PI), and Spring Tmax anomaly|Log(PI). We chose to test the two-way interactions individually rather than to include all of them in a single model in order to facilitate the biological interpretation of the coefficients of each interaction and main effect. Additionally, potential collinearity between predictors (such as the cross-products that comprise interactions) can lead to variance inflation, making detection of significant effects difficult (i.e., increased Type II errors).

Among all of these models, the only significant two-way interaction was the AP normal|Tmax normal interaction (see Results). This interaction term was also statistically significant (P = 0.0178) when tested in a model that included the six main effects above plus all of  two-, three-, and four-way interaction terms between and among the four climate variables (AP and Spring Tmax normal and anomalies; results not shown). We then tested for differences between species with respect to this interaction by constructing and testing a linear model that included the following predictors: Species, Log(PI), AP and Spring Tmax normals, AP and Spring Tmax anomalies, and the AP normal|Tmax normal and the Species|MAP normal|Tmax normal interaction terms.

Long-term phenological shifts. To detect long-term phenological shifts in flowering dates in each species across the sampling period, we constructed and tested linear models in which DOY was the response variable and log(PI), year, latitude, longitude, and elevation were treated as fixed independent variables. In these main-effects models, there was no significant effect of year on DOY. We then tested for two-way interactions between year and each geographic attribute: Year|Latitude, Year|Longitude, and Year|Elevation. Each two-way interaction was tested separately by adding it to the main-effects model. In the absence of any significant two-way interaction, a significant (or non-significant) effect of Year on DOY could be interpreted as an effect of Year that is independent of the effects on DOY of latitude, longitude, or elevation. In turn, a significant interaction would indicate that the rate of change in DOY over time (i.e., the effect of Year) varies among specimens collected in different latitudes, longitudes, or elevations.

All linear models were conducted using the lm function (stats 4.0.2 of the R Stats Package) in R Studio version 1.2.5042; figures were created using visreg v2.7.0 (Breheny and Burchett, 2017), and ggplot2 v3.3.2 (Wickham 2016). In all of these linear models, significance testing was conducted using Type III sums of squares; the effects on DOY of each independent variable or interaction term was tested when placed last into the model. In multiple linear models that included species as an independent variable, the lsmeans and eff_size functions (in the emmeans package) were used to test for significant differences between species’ least squares means and for effect sizes.

Please refer to the README.txt file.