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 yr; 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 d 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.
The relationship between phenological events — particularly the onset of flowering — and local climatic conditions has recently attracted the intense interest of ecologists and evolutionary biologists for several reasons (Wolkovich et al. 2012; Willis et al. 2017a). First, quantitative estimates of the sensitivity of flowering date to geographic variation in chronic climatic conditions or to inter-annual variation in climate over the last century can be used to forecast species-specific shifts in flowering date in response to upcoming climate warming. These predicted phenological shifts may, in turn, be used to infer changes in the exposure of flower buds or newly opened flowers to environmental stressors such as early spring frost (Park et al. 2020). Second, observed differences in sensitivity to climatic changes among sympatric or regionally co-occurring species lead to the expectation that community-level flowering patterns will change in response to a changing climate (Park and Mazer 2019). Similarly, differences in the sensitivities of co-occurring native and invasive plants provide insights into the causes of their relative success (Willis et al. 2010). Third, phenological shifts that result in asynchronous activity between mutualistic species (e.g., plants and their pollinators) or between plants and their antagonists (e.g., herbivores, florivores, or fungal spore vectors) may affect individual fitness or population mean fitness of all interacting taxa, and such mismatches may have cascading effects throughout an ecosystem (Kudo and Ida 2013; Kudo and Cooper 2019).
Detecting species-specific or multi-species phenological responses to climatic conditions requires observations that capture variation in both phenological behavior (comprising the response variable) and climatic parameters or other predictor variables, such as elevation, latitude, and longitude, that may serve as a proxy for climate. Given this requirement, herbarium collections are highly suited to the study of phenological variation and its environmental correlates for species that have been widely collected over space and time (Yost et al. 2019; Pearson 2019; Pearson et al. 2020, 2021). This suitability applies particularly well to annual forbs, for which herbarium specimens often contain intact individuals whose phenological status may be scored in a highly quantitative manner, thereby capturing the precise phenological condition of the entire individual plant when collected (Willis et al. 2017; Love et al. 2019; Goëau et al. 2020).
To date, most herbarium-based studies of phenological sensitivity to climate fall into one of three non-mutually exclusive categories. The first is represented by detailed studies of individual species in which the sensitivity to long-term mean climatic conditions and/or to inter-annual variation in climate is estimated, and some report interactions between predictor variables that affect phenological behavior (Robbirt et al. 2011; Gaira et al. 2011, 2014; Matthews and Mazer 2016; Ellwood et al. 2019; Love et al. 2019; Petrauski et al. 2019; Banaszak et al. 2020; Pearson et al. 2021). The second category comprises synthetic studies of multiple species and higher taxa, aiming to detect general similarities and differences among taxa or communities with respect to their phenological responses to climatic factors that vary over time or space (Primack et al. 2004; Miller-Rushing et al. 2006; Houle 2007; Gallagher et al. 2009; Diez et al. 2012; Panchen et al. 2012, 2017; Diskin et al. 2012; Calinger et al. 2013; Li et al. 2013; Mazer et al. 2013; Hart et al. 2014; Park 2014; Davis et al. 2015; Kharouba and Vellend 2015; Rawal et al. 2015; Munson and Long 2017; Park and Schwartz 2018; Jones and Daehler 2018; Park et al. 2018; Park and Mazer 2018, 2019; Berg et al. 2019; Pearson 2019; Kopp et al. 2020; Park, Ramirez Parada, and Mazer 2020; Reeb et al. 2020). A few studies in both categories have begun to investigate sources of intraspecific variation in phenological sensitivity (Matthews and Mazer 2016; Park et al. 2018; Song et al. 2020). The third category is represented by studies that compare the sensitivities derived from herbarium specimen data to those derived from high-resolution field observations (Robbirt et al. 2011; Davis et al. 2015; Zohner and Renner 2014) or from remotely sensed data (Park 2012) to determine the reliability of the former, which are known to include sampling biases of various kinds (taxonomic, geographic, seasonal, and temporal; Daru et al. 2017).
Few herbarium-based studies, however, have compared phenological sensitivities to climate of closely related species that are broadly sympatric and that share pollinators. Similarities between such species in phenological sensitivity to local climatic conditions might be expected due to several non-mutually exclusive processes or factors, including: (a) phylogenetic niche conservatism, whereby closely related species “inherit” their exposure to similar environmental conditions, resulting in parallel evolution of mean flowering date and/or sensitivity to the same environmental cues used to induce flowering; (b) the evolution of synchronous flowering because of its advantage in attracting a shared set of pollinators or in satiating predators; or (c) a lack of genetic variation in phenological sensitivity, precluding its evolutionary divergence between taxa. The use of herbarium specimens to compare congeners with respect to flowering time and climate-driven phenological sensitivity therefore provides an opportunity to address several fundamental ecological and evolutionary questions, including the following:
(1) Do congeners differ in the timing of flowering, controlling for spatiotemporal variation in local environmental variables that might influence flowering time, such as elevation, temperature, and precipitation? Closely related species that share pollinators and that exhibit overlapping geographic distributions may be expected to evolve to differ in mean flowering time when exposed to the same climatic conditions as a result of divergent natural selection favoring genotypes that avoid direct competition for pollinators or that minimize heterospecific pollen transfer (Rathcke and Lacey 1985; Fenner 1998). In such cases, divergent flowering times would contribute to reproductive isolation, reducing the potential for hybridization and/or heterospecific pollination, both of which may reduce seed production or quality (Morales and Traveset 2008; Moreira-Hernández and Muchhala 2019; Ashman and Arceo-Gómez 2013; Arceo-Gómez et al. 2019), particularly when heterospecific pollen transfer occurs between close relatives (Streher et al. 2020).
(2) Independent of mean flowering times, do closely related taxa exhibit distinct sensitivities of flowering time to local temperature and precipitation? Close relatives that share habitats and pollinators might be expected to evolve contrasting or opposing phenological responsiveness to a given set of environmental cues because, just as is the case for mean flowering date, divergence in sensitivity may promote reproductive isolation and reduce heterospecific pollen transfer. Alternatively, species that share pollinators might be expected to evolve to use the same environmental cues to induce flowering where synchronous flowering attracts more pollinators per capita or per flower due to potential synergistic effects of producing a large floral display at the community level (McGuire 1993; Kudo 2006; Sargent and Ackerly 2008).
(3) Does the flowering phenology of closely related species respond similarly to geographic variation in long-term mean climatic conditions and to inter-annual variation in local climate? The process that generates associations between flowering date and long-term mean climatic conditions may differ from the process (or processes) that causes flowering date to respond to inter-annual variation in climate. Strong associations between long-term, chronic, local climatic conditions and flowering date likely represent the outcome of natural selection, while phenological shifts in response to comparatively unpredictable inter-annual variation may represent adaptive or non-adaptive plastic responses to short-term environmental variation or, potentially, the effects of rapid adaptive evolution. Given that different processes (adaptation vs. plasticity) may contribute to phenological sensitivity to long-term vs. inter-annual climatic conditions, both the magnitude and the direction of these sensitivities may differ within or between species (Mazer et al. 2020).
In the current study, we address these questions using physical herbarium specimens of Clarkia cylindrica and C. unguiculata (Onagraceae), a pair of ecologically similar annual herbs native to California. To determine whether these Clarkia congeners respond similarly to climate or to climate change, we tested for differences between them with respect to the climatic conditions they occupied; estimated the degree of climate change they experienced in their sampled ranges during the past century; used phenoclimatic models to estimate the phenological responsiveness (i.e., sensitivity) of each species to both long-term mean climate and interannual variation in climate; and tested for long-term shifts in the estimated flowering date over the past century in both species.
Methods
Overview
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 yr across its 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 (2021) investigation of Eschscholzia californica (California Poppy; Papaveraceae), and Parker's (2021) 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 normals 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-yr 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 a 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 (section Peripetasma) and C. unguiculata Lindley (section 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 of the latter 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 1500 m, 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 (>5 mm 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; Wang et al. 2016). 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 <707 m from another collection site. For 85 specimens of Clarkia cylindrica, the distance to the nearest collection site was <707 m, 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 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 <707 m 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. 2019): 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 × 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 R2 values) than those that included Spring Tmin ( Appendix S1 (388_madr-68-04-08_s01.docx)). 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 yr 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-yr 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 × Spring Tmax normal, AP normal × AP anomaly, AP normal × Spring Tmax anomaly, Spring Tmax normal × AP anomaly, Spring Tmax normal × Spring Tmax anomaly, AP anomaly × Spring 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) when models contain all possible predictors.
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 the 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) differs 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.
Results
Climatic Conditions and Geographic Locations Occupied by Each Species
Clarkia cylindrica and C. unguiculata are endemic to California, with C. unguiculata being the more geographically widespread of the two (Fig. 1). Bivariate plots illustrating the joint distributions of the sampled specimens of the two species with respect AP normal, spring Tmax normals, elevation, longitude and latitude reveal broad overlap between the two species with respect to all variables, although C. unguiculata extends farther to the north, west, and southeast of California than C. cylindrica ( Fig. S1 (388_madr-68-04-08_s01.docx)).
The linear models designed to test for differences between the two species with respect to their exposure to long-term variation in climate detected that the specimens of C. cylindrica were collected from chronically drier and cooler sites than those of C. unguiculata (Fig. 2). The least squares mean (LSM) for the AP normal from 1921–2010 was 463 mm (±SE = 10.3) for C. cylindrica and 563 mm (±SE = 11.5) for C. unguiculata (effect size = 0.58 ± SE = 0.09, df = 504, P < 0.0001). The LSM for the Spring Tmax normal was 19.8°C (±SE = 0.107) for C. cylindrica and 20.4°C (±SE = 0.120) for C. unguiculata (effect size = 0.321 ± 0.92, df = 504, P = 0.0005). The collection sites of the two species, however, exhibit similar and near-zero mean values of the AP and spring Tmax anomalies (Fig. 2), indicating that neither species was disproportionately represented by sites that were hotter (or cooler) or drier (or wetter) than average in the year of collection.
Fig. 2.
Box plots displaying the least squares means (LSMs) for C. cylindrica and C. unguiculata for each of four climatic variables, where the LSM of the focal dependent variable for each species was estimated while controlling for variation in the other three variables. Based on the specimens sampled here, the two species differ significantly with respect to the mean values of the chronic climatic conditions (mean annual values from 1921–2010) for cumulative annual precipitation (AP) and Spring Tmax. C. cylindrica occupies drier and cooler sites than C. unguiculata. Horizontal black lines represent the median of each climatic variable; diamonds represent the mean values. Within each panel, if the outer openings of the notches of the two species do not overlap, this can be interpreted as a significant difference between the median values of the two species.
The mean elevation of the sites from which C. cylindrica was collected was significantly higher than that of the sites from which C. unguiculata was sampled (one-way ANOVA of log-transformed values of elevation: F1, 507 = 25.95, P < 0.0001, R2 = 0.05). Clarkia cylindrica specimens were sampled at a mean elevation of 636.3 m (±SE = 21.90, n = 283) and C. unguiculata specimens were collected at a mean elevation of 493.4 m (±SE = 23.4, n = 226) ( Fig. S2 (388_madr-68-04-08_s02.docx)).
Climate Change During Sampling Period: Temporal Change in Climate Anomalies
Both species were exposed to similar patterns of temporal change in climate anomalies over the ∼115-yr sampling periods with respect to AP anomalies, winter and spring Tmax and Tmin anomalies, and the anomalies for the number of frost-free days (NFFD) in winter and spring. Both species experienced small, but statistically significant, increases in AP anomalies (90.4 and 114 mm/100 yrs, respectively for C. cylindrica and C. unguiculata) (Table 1). Similarly, both species experienced significant temporal increases in winter and spring Tmin anomalies (Fig. 3), but no significant change in either winter or spring Tmax anomalies over the sampling period. Finally, both species experienced similar and statistically significant temporal increases in the number of frost-free days in both winter and spring (Fig. 3).
Table 1.
Summary of Bivariate Regressions Conducted Using the Specimens of Each Species Separately to Detect Temporal Trends in Climate Anomalies Over the ∼115-year Sampling Period. The two species show similar magnitudes of temporal change in the anomalies for cumulative annual precipitation (AP, including rain and snow), PAS (precipitation as snow), winter and spring mean minimum daily temperature (Tmin), and winter and spring NFFD (number of frost-free days), indicating small, but significant increases in precipitation and winter and spring temperatures. Boldfaced values indicate variables for which significant change (α = 0.05) was detected over the sampling period (values are italicized when 0.05 < P < 0.10). The inclusion of elevation, latitude, and longitude as main effects in these models did not significantly change the magnitude of the parameter estimate for the effect of year on any climate variable.
Fig. 3.
Bivariate regressions illustrating significant temporal trends in winter and spring Tmin anomalies (the deviations between the climate normals and conditions in the year of specimen collection at each site) and in the anomalies for the number of frost-free days (NNFD) in winter and spring. See Table 1 for regression parameters.
Factors Influencing Flowering Date: Species, Spring Tmax, and AP Normals and Anomalies, and Phenological Status
Differences between species' mean DOY. Controlling for the effects on DOY of the other predictor variables in the model, C. cylindrica flowers 3.46 d earlier (±SE = 1.55; effect size =–0.21 ± SE = 0.932, df = 501) than C. unguiculata (LSM of DOY for C. cylindrica = 146.0 ± SE = 1.03; LSM for C. unguiculata = 149.0 ± SE = 1.15; t-ratio = 2.232, P = 0.0260) (Table 2; Fig. S3 (388_madr-68-04-08_s03.docx)).
Table 2.
Summary of Mean-Centered Linear Model Designed to Detect Significant Effects on the Day of Year (DOY) of Specimen Collection of Phenological Status, Species, Geographic and Inter-Annual Variation in Climate, and the Two-Way Interaction Between Long-Term Annual Precipitation (AP Normal) and Long-Term Spring Tmax (Mean Maximum Daily Temperature). The values of the phenological index were log-transformed to improve normality. No significant interactions were detected between Species × AP normal, Species × Spring Tmax normal, AP normal × Spring Tmax anomaly, or Species × AP normal × Spring Tmax normal, so these interactions are excluded here. The effects on DOY of the main effects and interactions shown here do not differ between species. Type III sums of squares were used for significance testing.
Effect of long-term, chronic climatic conditions. The linear model applied to the data set in which the two species were pooled detected significant effects of AP normals and Spring Tmax normals on DOY while controlling for the phenological status of specimens, inter-annual variation in AP and spring Tmax, and the two-way interaction between AP normals and spring Tmax normals (Table 2).
The significant two-way interaction between AP and Spring Tmax normals indicates that the delaying effects of high precipitation are stronger in warmer than in cooler regions and that the advancing effects of higher spring Tmax are stronger in chronically dry than in relatively wet conditions (Fig. 4). Flowering consistently occurs earliest under hot and dry conditions. This interaction did not differ significantly between species ( Appendix S2 (388_madr-68-04-08_s02.docx); the Species × AP normal × Spring Tmax normal interaction was non-significant) although when we tested species-specific models, a significant interaction was detected in C. cylindrica but not in C. unguiculata ( Fig. S4 (388_madr-68-04-08_s04.docx), Appendix S3 (388_madr-68-04-08_s03.docx)). We detected no evidence in the analysis of the pooled data set that species differed in their sensitivity to AP or Spring Tmax normals. When added to the model containing only main effects (Species, PI, AP normals and anomalies, and Spring Tmax normals and anomalies), no significant two-way interactions were detected between Species and either AP or Spring Tmax normals (results not shown).
Fig. 4.
Illustration of the significant interaction between long-term mean annual precipitation (Annual Precipitation Normal) and long-term mean maximum spring temperature (Spring Tmax normal) detected in the analysis of pooled data (Table 2). A. The delaying effects on DOY of high precipitation are strongest in relatively warm regions. B. The advancing effects on DOY of higher mean Spring Tmax are strongest in relatively dry regions.
Effects of inter-annual variation in AP and Spring Tmax. Among all specimens, we detected a 1.0-d delay for every 100 mm increase in the AP anomaly and a 4.2-d advancement for every 1°C increase in the Spring Tmax anomaly. We detected no significant two-way interaction between the AP normal and Spring Tmax anomalies (Table 2). Moreover, we detected no evidence in the analysis of the pooled data set that species differed in their sensitivity to AP or Spring Tmax anomalies. When added to the model containing only main effects, no significant two-way interactions were detected between Species and either AP or Spring Tmax anomalies (results not shown).
Effects of phenological status on flowering date. When added to the model containing only main effects, we detected no significant PI × Species interaction, indicating that the rate of phenological progression over time is the same for both species (Fig. 5). The necessarily positive slope of the relationship between DOY and log(PI) (i.e., the regression coefficient associated with log[PI] in Table 2) indicates how rapidly plants progress from bearing all buds to bearing only fruits. Higher values of this slope indicate that it takes longer to progress between phenological stages; lower values (flatter lines) indicate that progression between stages occurs more rapidly.
Fig. 5.
Leverage plots illustrating the relationship between day of year of specimen collection and the phenological status of herbarium specimens representing Clarkia cylindrica and C. unguiculata. Herbarium specimens represented by a value of 0 for the Log(Phenological Index) bear only buds; those represented by a value of 0.6 bear only ripe fruits. The absence of a significant Species × Log(PI) interaction ( Table S5 (388_madr-68-04-08_s05.docx)) indicates that the slope of the relationship shown here (67.29 ± SE = 6.42) does not differ between species.
Just as the rate of phenological progression did not differ between species, the rate of progression was also unaffected by the climatic variables examined here. None of the following two-way interactions were significant: PI × AP normal, PI × Spring Tmax normal, PI × AP anomaly, or PI × Spring Tmax anomaly (results not shown), indicating that the rate of phenological progression did not depend on local climatic conditions.
Long-term Phenological Shifts
Neither species exhibited statistically significant phenological change over the ∼115-yr sampling period, independent of the effects on DOY of phenological status, elevation, latitude and longitude (Table 3). These models detected, however, that in both species DOY changes significantly with elevation and longitude; higher elevations and more western locations are associated with delayed flowering (Fig. 6). The species-specific models containing only main effects (Table 3) have the same adjusted model R2 values (0.23 for C. cylindrica; 0.19 for C. unguiculata) as the three species-specific models that each contained a two-way interaction (Year × Elevation, Year × Latitude, or Year × Longitude), none of which were significant ( Table S4 (388_madr-68-04-08_s04.docx): Adjusted R2 values for models including the Year × Elevation, Year × Latitude, and Year × Longitude interaction term, respectively, are 0.23, 0.23, and 0.24 for C. cylindrica, and 0.19, 0.19 and 0.19 for C. unguiculata).
Table 3.
Summary of Species-Specific Parameter Estimates (Slopes) for Mean-Centered Linear Models Designed to Detect Long-Term Temporal Trends in Flowering Date (DOY) While Controlling for Variation in Phenological Status (Log[PI]), Elevation, Latitude, and Longitude. A. Clarkia cylindrica. B. C. unguiculata. Two-way interactions involving year (Year × Elevation, Year × Latitude, and Year × Longitude) were all non-significant (see Table S4 (388_madr-68-04-08_s04.docx)), supporting the interpretation that the effects of year on the day of year of specimen collection are consistent across elevation, latitude, and longitude. Type III sums of squares were used for significance testing.
Fig. 6.
Leverage plots illustrating the relationship between day of year of specimen collection and elevation (A and C) and longitude (B and D), controlling for variation in the phenological status of specimens, year of collection, and latitude. See Table 3 for parameter estimates (slopes, standard errors, and P-values).
Discussion
Among the specimens sampled here, C. cylindrica and C. unguiculata exhibit more similarities than differences in their phenological behavior. They do not differ in their sensitivities to chronic local climatic conditions (normals) or to year-specific conditions (anomalies). In addition, they exhibit the same rates of progression through successive phenological stages (Fig. 5), which appear to be constant throughout their ranges (there were no significant interactions between PI and any of the four focal climate variables). If these similarities are in fact the result of adaptive evolution, then the benefits of co-flowering would appear to have outweighed the benefits of divergence and reproductive isolation. These similarities predict that, all else being equal, as both precipitation and spring temperature change in the future, the two species' flowering phenology should shift to the same degree and in the same direction.
Difference Between Species in the Timing of Flowering
The specimens of C. cylindrica and C. unguiculata sampled here differ with respect to the historical conditions at the sites from which they were collected, with C. cylindrica occupying drier and cooler locations than C. unguiculata (Fig. 2). This difference between the species in the mean climatic conditions of the sampled sites from 1920–2010, however, does not explain the statistically significant 3.5-d difference between them in DOY, because independent of local climatic conditions (both the normals and the anomalies) and the phenological status of specimens, this 3.5-d difference persists (Table 2). Whether this difference between species in mean DOY represents the adaptive outcome of natural selection due to species-specific effects of flowering time on pollination, herbivory, or the intensity of inter-specific competition, or whether it represents an evolutionary or plastic response to other environmental differences between the species that were not investigated here, cannot be determined without additional research. More fundamentally, it is not clear that a difference of this magnitude has any ecological or fitness consequences.
Species Exhibit Similar Sensitivities of Flowering Time to Local Temperature and Precipitation, and Similar Rates of Phenological Progression
The absence of statistically significant interactions between species and the focal climate variables (Species × AP normal, Species × Spring Tmax normal, Species × AP anomaly, Species × Spring Tmax anomaly, and Species × MAP normal × Spring Tmax normal; results not shown) indicates that species do not differ with respect to the direction or magnitude of their responses to local AP or Spring Tmax normals or to warmer- or wetter-than average conditions in the year of collection (represented by the AP and Spring Tmax anomalies). In addition, the lack of a significant Species × Log(PI) interaction ( Table S5 (388_madr-68-04-08_s05.docx)) and the parallel leverage plots of the relationship between DOY and log(PI) (Fig. 5) indicate that the rates at which the two species progress through their successive phenological stages do not differ, based on these analyses.
Other studies of congeners have observed both similarities and differences in their phenological behavior, in part depending on whether the species' native ranges encompassed the same climatic conditions. In a study of four Quercus species — two species restricted to a water-limited, Mediterranean climate in California (Q. agrifolia Née and Q. lobata Née) and two species occupying the temperate zone of eastern and central North America (Q. alba L. and Q. rubra L.) — Gerst et al. (2017) found that individual trees of the western species tended to flower intermittently throughout the growing season and for much longer durations than the eastern/central species, which exhibited much more synchronous and fewer flowering onset dates. Gerst et al. (2017) interpreted the western species' behavior as a potential adaptive evolutionary or plastic response to unpredictable episodes of rainfall throughout the late winter and spring in California. The California oaks appeared to flower opportunistically in response to these pulses of precipitation, while the eastern/central U.S. species tracked the transition between winter and spring, and flowered comparatively synchronously in response to local winter or spring Tmax.
A similar pattern was found among 11 congeneric pairs of high- vs. low-elevation species of perennial herbs in Switzerland; the elevation of origin determined the sensitivity of phenological traits to local growing conditions (Schmid et al. 2017). Each of the 11 species pairs were grown in both a high- and a low-elevation common garden and a suite of phenological events was recorded. The species native to low-elevations showed much stronger plastic responses of flowering phenology to the growing conditions of the common gardens than their high-elevation counterparts, suggesting that phenological sensitivity is a trait that may evolve in response to local conditions.
An experimental study of two temperate zone congeners of Echium (Boraginaceae) detected similar phenological sensitivity to drought and heat stress between the short-lived perennial, E. plantagineum L., a native of Southern Europe and the annual, E. vulgare L., native to Northern Europe (Descamps et al. 2020). In both taxa, experimentally induced stress advanced flowering time, but the rate at which the two species produced sequential flowers differed, with the perennial species exhibiting stronger phenological responses to the distinct watering and heat treatments than the annual.
Perhaps the most robust way to seek evidence for evolutionary similarities due to the benefits of sharing pollinators (facilitation) vs. the benefits of phenological divergence (reproductive isolation) is to compare the flowering times of congeners when they appear in allopatry vs. sympatry. In a comparative study of 74 pairs of congeners in 26 genera, Park et al. (2020) found that flowering time displacement between congeners observed in sympatry was uncommon, consistent with the hypothesis that the benefits of reproductive isolation are relatively small. Several additional broad, comparative studies have also demonstrated that closely related species tend to have similar flowering dates or sensitivities (Davies et al. 2013 [∼4000 species in 1246 genera]; Du et al. 2015 [19631 species in 2284 genera and 181 families], 2020 [2136 species in 846 genera and 158 families]). These studies corroborate the hypothesis that phylogenetic niche conservatism (Wiens and Graham 2005), genetically based trait conservatism (Harvey and Pagel 1991), and parallel or convergent evolution of flowering date (where synchronous flowering facilitates pollinator-sharing) may be more common than selection favoring phenological divergence in flowering time between close relatives. Considered together, these studies suggest that the benefits of co-flowering for pollinator attraction and visitation are greater than the benefits of reproductive isolation.
Sensitivity to Long-Term Climatic Conditions vs. to Inter-Annual Variation in Climate
Higher precipitation and cooler temperatures delay flowering in C. cylindrica and C. unguiculata, whether those climatic conditions vary geographically or interannually (Table 2). The two species exhibit similar sensitivities to long-term AP normals and AP anomalies (0.018 ± SE = 0.005 vs. 0.010 ± 0.004 d/mm, respectively) and to long-term spring Tmax normals and anomalies (–2.75 ± 0.42 vs. –4.19 ± 0.60 d/°C, respectively). For these species and climatic variables, therefore, estimates of sensitivity to spatial variation in climate appear to provide reasonable “space-for-time” proxies that can be used to forecast the direction of long-term phenological shifts in response to upcoming change in these parameters. Nevertheless, in the data analyzed here, these sensitivities were not consistent with the absence of statistically significant phenological shifts observed over time.
The linear models tested here include as predictor variables both the climate normals (variation among which represents geographic variation in mean climate conditions) and the climate anomalies (which result from interannual variation in climate) for each site of specimen collection. Accordingly, the effects on DOY of each climatic variable are independent of the effects of the others (Table 2). The significant effects on DOY detected for both the normals and the anomalies are therefore consistent with the interpretation that both evolutionary adaptation to local temperature and precipitation, and plasticity induced by these parameters in the year of collection, contribute to variation in DOY in these species. In addition, we cannot rule out the possibility that recent, short-term evolutionary responses to temperature and precipitation anomalies also contribute to the associations between DOY and these anomalies (Table 2). Papper and Ackerly (2021) similarly found evidence for both local adaptation and phenotypic plasticity of the phenology of bud break in Quercus douglasii Hook. & Arn. (Fagaceae), demonstrating that the ability to distinguish between these processes is not restricted to short-lived annual species.
Current end-of-century climate projections in California predict statewide increases in temperature that range from 2–4°C (assuming an intermediate level of emissions) and from 4–7°C (assuming no mitigation of emissions), while predicted changes in annual precipitation are regional, season-specific, and less certain, as inter-annual variation in precipitation is expected to increase (Pierce et al. 2018). Assuming that future phenological shifts in these taxa can be predicted from observed responses to inter-annual variation in mean spring maximum temperature (i.e., to the anomalies) alone, we would predict that upcoming warming in California will generate phenological advances in their flowering that range from ∼8.4–28 d by the end of the century. The phenological response to joint changes in both temperature and precipitation, however, as well as the consequences of these shifts for altering the exposure of wild Clarkia populations to herbivores, to pathogens, and to pollinators are obviously unknown.
Long-term Phenological Shifts
Although both species exhibit significant phenological advancement in response to warmer-than-average years (–4.19 d/°C increase in spring Tmax; Table 2), this sensitivity did not generate significant phenological change in either species over the ∼115-yr sampling period (Table 3; the effect of year on DOY was not significant). The simplest explanation for this is that there was not a sufficient increase in Spring Tmax over the sampling period in either species to induce detectable shifts in DOY over time (Table 1). However, Spring Tmin anomalies did increase significantly over time (1.27°C/century for C. cylindrica and 1.45°C/century for C. unguiculata: Table 1), and DOY was sensitive to increases in Spring Tmin anomalies (advancing 3.2 d/°C; Table A6 (388_madr-68-04-08_s06.docx)). These values predict an advancement in DOY of 4.1 d/century in C. cylindrica and 4.6 d/century in C. unguiculata, which were not observed. Other factors that affect DOY and that changed over the sampling period but were not examined here may have opposed this predicted advancement in flowering.
One factor that could have opposed a long-term advancement in DOY is the temporal increase in annual precipitation anomalies (90.4 mm/century for C. cylindrica and 114 mm/century for C. unguiculata: Table 1). Given the sensitivity of these species to AP anomalies (a delay in DOY of 0.01 d/mm increase in the anomaly; Table 2), the temporal increase in AP over the sampling period is expected to have generated a 0.90-d/century delay in flowering in C. cylindrica and a 1.14-d/century delay in flowering in C. unguiculata. These predicted delays, however, would not have been sufficient to offset the predicted advancement due to warmer spring temperatures.
In a study of five chaparral species of shrubs in the genera Arctostaphylos (Ericaeae) and Ceanothus (Rhamnaceae), Parker (2021) similarly detected no significant temporal shift in the flowering phenology even though the mean annual temperature of the collection sites increased significantly across the collection period. By contrast, Pearson et al. (2021) detected a significant advancement in flowering time of Eschscholzia californica of 0.6 d per decade across its range, with the magnitude of advancement varying with local long-term mean annual temperature.
The models testing for long-term phenological shifts independent of variation in elevation, latitude, and longitude (Table 3) detect that in both species DOY changes significantly with elevation and longitude; higher elevations and more western locations are associated with delayed flowering (Fig. 6). The non-significant interaction terms between year and each of the geographic variables ( Table S4 (388_madr-68-04-08_s04.docx)) indicate that the temporal trends in phenology did not depend on elevation, latitude or longitude. Therefore, the non-significant effect of year cannot be attributed to biased sampling of elevations or longitudes over time. In other words, higher elevations or more western sites might have been disproportionately sampled in more recent years, cancelling out the advancing effects of warming on DOY, but this is not the case.
Selection on Flowering Date in Clarkia: Empirical Evidence from Other Studies
Several studies of wild populations of Clarkia species corroborate the view that sympatric populations benefit from the presence of congeners by attracting pollinators and facilitating each other's pollination. Moeller (2005) observed that, in the southern California endemic, C. xantiana A.Gray subspecies xantiana, populations growing in sympatry with congeners (including C. cylindrica and C. unguiculata) with which they share pollinators (the most effective of which he reported to be Hesperapis regularis, Lasioglossum pullilabre, and Megachile gravita or M. pascoensis) had visitation rates that were two times higher than those of populations growing without congeners (see also Moeller and Geber 2005). In addition, congeneric populations growing in sympatry had lower degrees of pollen-limited reproduction than those growing in isolation (Moeller 2004). Moreover, populations of C. xantiana subsp. xantiana growing in sympatry with other congeners did not differ in their mean flowering time from those growing alone, indicating that there is no evidence of character displacement in this trait in response to inter-specific interactions (including possible competition for pollinators) (Moeller 2004).
One may, with caution, use herbarium specimens representing well-sampled taxa to detect the evolutionary effects on flowering time of growing in sympatry or allopatry with congeners (Park et al. 2020). In the current study, due to both the broadly overlapping ranges of C. unguiculata and C. cylindrica, and the small number of specimens of each taxon collected from sites that appear to be beyond the range of the other taxon (Fig. 1), it would be risky to use the data analyzed here to reliably identify sites of long-term allopatry vs. sympatry. Moreover, specimen-based data indicate only whether a given species was present at a given location and time, but not if it was absent, making it difficult to assert allopatry in the absence of direct field observations. Nevertheless, herbarium-derived data for well-sampled taxa could be used to identify potential areas of allopatry or sympatry that could then be surveyed in the field to confirm presence or absence of each species. Once these areas are identified, phenological data derived from herbarium specimens could be used to address whether, independent of climate, the mean phenological difference between taxa depends on whether they are in sympatry vs. allopatry.
Other evidence supports the view that competition for pollinators does not select for phenological divergence between sympatric Clarkia congeners. Eisen et al. (2019a) conducted a field experiment in which the flowering times and the number of co-occurring Clarkia species were manipulated by placing potted, late-flowering plants in natural communities with or without earlier-flowering congeners. They found that seed set of the relatively late-flowering taxon, C. xantiana subsp. xantiana, depended on local floral density, but that its delayed flowering relative to other congeners could not be attributed to selection favoring reduced interspecific competition for pollinators. Eisen et al. (2019b) also report that the relatively late flowering of C. xantiana subsp. xantiana is not adaptive; plants had higher seed set when they flowered earlier. They proposed, therefore, that the staggered flowering times among co-occurring Clarkia taxa may be due to ecological specialization unrelated to pollination, to developmental constraints associated with other components of their life history (e.g., the need to reach a threshold size before flowering), or to genetically based correlations with other traits that are under strong selection and drive the evolution of species-specific flowering times as a consequence.
If these patterns apply to Clarkia unguiculata and C. cylindrica, then the phenological similarities between them reported here may reflect phylogenetic niche conservatism (cf. Truszczynski et al. 2021), trait conservatism, or the outcome of selection on flowering time due to facilitation between co-flowering species that share pollinators. It must be emphasized, however, that although herbarium specimen-based studies can detect phenological patterns over a broader geographic and temporal range than short-term field studies, they are limited in their capacity to identify the mechanisms that contribute to observed phenological divergence or convergence, underscoring the need for complementary studies of natural populations. To detect the factors that independently and directly promote synchronous flowering between congeners, additional observational and experimental work in field populations and communities are needed to measure the costs and benefits associated with synchronous vs. asynchronous flowering between congeners. The continued analysis of large herbarium-based data sets that include congeners, however, hold the promise of detecting those phenological similarities and differences between congeners that merit investigation in the first place.
Data Accessibility
The data sets and R code used in this paper are archived in the Dryad Digital Repository ( https://doi.org/10.25349/D9N026).
Acknowledgments
The authors gratefully acknowledge the many UCSB students who invested much time and effort in the scoring and georeferencing of hundreds of specimens, including Adam Lota, Adam Monahan, Allison White, Dana Hoffenberg, Elizabeth Cannon, Emmanuel Lopez, Jazmine Brown-lee, Joel Kirksey, Kimberly Crispin, Louis Lin, Laurel Schmalz, Marco Marquez, Melissa Mckiddy, Robert Ottolia, Tahirany Diaz, and Tiffany Thaphasouk. We also thank Danica Taber for her assistance in managing students and data. We are very grateful to the staff of the JEPS, RAS, SBBG, UC, and UCR herbaria who managed the loans for this project, and particularly to Dr. Jennifer Thorsch and the herbarium of UCSB's Cheadle Center for Biodiversity and Ecological Restoration for hosting the borrowed specimens and for providing workspace. The authors acknowledge funding from NSF (DBI-1802181) to SJM, from the National Park Service Climate Change Response Program (Cooperative Agreement H8C07080001), and from UCSB, which provided a generous Chancellor's Fellowship to TRP. SJM is indebted to the Yale Institute for Biospheric Studies, which provided support and camaraderie during the sabbatical year in which this manuscript was prepared. Finally, this work would not have been possible without the generosity and time of the many collectors and systematists who collected and identified the specimens analyzed here.
