DOI: 10.25349/d9n026
title: Data used in: Phenological sensitivities to climate are similar in two Clarkia congeners: Indirect evidence for facilitation, convergence, niche conservatism, or genetic constraints

description: 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.
description: 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 <i>C.
 cylindrica</i> and 1892-2011 for <i>C. unguiculata</i>).
 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
 <i>deviation</i> between climatic conditions at the site in
 the year of collection and the site’s mean long-term climate conditions
 (i.e., climate <i>anomalies</i> 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 <i>Clarkia</i>; Pearson et al.’s (2022)
 investigation of <i>Eschscholzia californica</i> (California
 Poppy; Papaveraceae), and Parker’s (2022) study of five species of
 <i>Arctostaphylos</i> (Ericaceae) and
 <i>Ceanothus</i> (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.
 <b>Study species</b>
 <i>Clarkia</i> 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
 <i>Clarkia</i> 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 –
 <i>Clarkia cylindrica</i> ssp. <i>clavicarpa</i>
 W. Davis (subgenus Peripetasma) and <i>C. unguiculata</i>
 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.
 <i>Clarkia cylindrica</i> has been described as “normally
 outcrossing” (Davis, 1970), and <i>C. unguiculata</i> 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 <i>Clarkia</i>, <i>C.
 unguiculata</i> and <i>C. cylindrica</i> could be found
 in mixed colonies with, respectively, any of 17 or 10 other
 <i>Clarkia</i> 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 <i>C.
 unguiculata</i> and <i>C. cylindrica</i> 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
 <i>C. cylindrica</i> and the rotate, clawed-petalled flowers
 of <i>C. unguiculata,</i> both of which provide pollen and
 nectar rewards (Fig. 1). <i>Clarkia cylindrica</i> is commonly
 visited for its pollen by <i>Andrena lewisorum</i> Thorp,
 which has also been found to visit <i>C. unguiculata</i>
 (MacSwain et al., 1973). Other bee species that visit <i>C.
 cylindrica</i> and <i>C. unguiculata</i> (but not
 necessarily at the same locations, at similar frequencies per plant
 species, or with equivalent pollination efficacy) include:
 <i>Andrena omninigra</i> Viereck, <i>Hesperapis
 regularis</i> (Cresson), <i>Megachile gravita</i>
 Mitchell, <i>M. pascoensis</i> Mitchell, <i>Diadasia
 angusticeps</i> Timberlake, <i>Melissodes clarkiae</i>
 LaBerge, <i>Synhalonia venusta carinata</i> (Timberlake),
 <i>Ceratina sequoiae</i> Michener, <i>Lasioglossum
 pullilabre</i> (Vachal) and <i>Apis mellifera</i>
 Linnaeus (MacSwain et al., 1973). <b>Herbarium
 data</b> The two <i>Clarkia</i>
 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 (&gt; 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 <i>Px</i> represents the proportion of reproductive organs in a given stage and <i>i</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 <i>C. unguiculata</i>, a total of 608 plants on 231 sheets were scored; for <i>C. cylindrica</i>, a total of 1042 plants on 306 sheets were scored. Duplicate specimens were considered to be those that were: collected &lt; 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 <i>C. unguiculata</i> collected from 1902-2011 and 284 sheets of <i>C. cylindrica</i> collected from 1900-2012 were analyzed here. Within each species, some specimens were retained even if they were collected &lt; 707m from another collection site. For 85 specimens of <i>Clarkia cylindrica</i>, the distance to the nearest collection site was &lt; 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 <i>Clarkia unguiculata</i> were collected &lt; 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 <i>C. cylindrica</i> was 1.0-3.77; the range for <i>C. unguiculata</i> was 1.0 – 3.95. In all analyses described below, the PI was log-transformed to improve normality.   <b>Climate data</b> 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 <i>x</i> 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 <i>C. cylindrica</i> and <i>C. unguiculata</i>, 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<sup>2</sup> 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 &lt; 0.0001, n=509; Spring Tmax anomaly vs. spring Tmin anomaly: r = 0.75, P &lt; 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.   <b>Analyses</b> <i>Climatic conditions and geographic locations occupied by each species. </i>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. <i>Climate change during sampling period: Temporal change in climate anomalies. </i>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. <i>Factors influencing flowering date: species identity, spring Tmax and AP normals and anomalies, and phenological status. </i>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. <i>Long-term phenological shifts. </i>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 <i>lm</i> function (stats 4.0.2 of the R Stats Package) in R Studio version 1.2.5042; figures were created using <i>visreg</i> v2.7.0 (Breheny and Burchett, 2017), and <i>ggplot2 </i>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 <i>lsmeans</i> and <i>eff_size</i> functions (in the <i>emmeans</i> package) were used to test for significant differences between species’ least squares means and for effect sizes.
description: Please refer to the README.txt file.