Concordance of long-term shifts with climate warming varies among phenological events and herbaceous species

. Many temperate herbs now flower earlier than a few decades ago. Little is known about other phenological events, despite the importance of life history integration for plant fitness. This study addresses the hypothesis that temporal shifts of multiple phenological events in herbs are concordant with temporal changes in weather. Explicitly showing that changes in timing of annual life cycle events are correlated with changes in weather-predicting variables provides support for the hypothesis that a phenological shift is concordant with climate change. Observations of six phenological events and five phenophases were made year- round for 25 yr for herb species in a deciduous forest fragment, Trelease Woods in Illinois, USA. Dates for 43 species were analyzed by linear mixed-effects models for events and pheno- phases and were compared to weather data from a nearby station. For early species, emergence was delayed by 1.5 d/decade, while end expansion advanced by 3.8 d/decade and begin dor- mancy advanced by 2.5 d/decade. For late species, end expansion advanced by 6.7 d/decade, while begin senescence delayed by 17.7 d/decade. Begin flowering and end flowering advanced similarly for both seasonal groups, at 3.8 to 4.2 d/decade. Some events showed no temporal change. Species differed greatly in the degree or direction of change, related to seasonality of event or length of phase. Overall, for a given species, most events are advancing (68.4%) and most durations are shortening (74.4%). In 12 of 13 cases, inter-annual variation in event date was predicted by a temperature-event – season combination variable, but in only six cases did both event date shift and weather variable warm through time. This finding supports the hypothesis that climate change is associated with changes in some, but not all, phenological events. This first long-term, multi-phase study of a community of temperate herb species indi- cates little temporal coherence of responses of multiple phases. Changes in date are event specific, phase specific, and species specific. This complexity of responses among species and uneven responses within a species ’ integrated annual cycle events has implications for evolutionary responses and more immediate interactions among plant, animal, and microbe species in this community.


INTRODUCTION
Long-term advances in dates of temperate phenological events have been linked to climate change, specifically rising temperature (Walther et al. 2002, Parmeson andYohe 2003). For temperate deciduous trees, climateassociated changes in spring leaf unfolding, autumn leaf coloration and drop, and flowering are well studied at local (Richardson et al. 2006), regional (Xie et al. 2015), and continental scales (Menzel 2000). In contrast, studies of herbaceous plant phenology are uncommon (Cleland et al. 2006) and multiple phases are rarely studied (Wolkovich and Ettinger 2014, Li et al. 2016, Heberling et al. 2019a. The focus of all long-term studies of the phenology of herbaceous communities has been on only flowering (Abu-Asab et al. 2001, Fitter and Fitter 2002, Inouye 2008, Calinger et al. 2013, in part, because many studies rely on historical records of early naturalists who restricted observations to a limited number of events (Bradley et al. 1999, Cook et al. 2008, Amano et al. 2010, Beaubien and Hamann 2011, Bolmgren et al. 2013, McDonough MacKenzie et al. 2019. Lacking are long-term, community-level, observational studies of herb phenology that incorporate all annual life cycle events. Short-term studies include species representing multiple syndromes with a focus on the link of phenology to carbon gain (Heberling et al. 2019a). Long-term phenological changes in the herb community have important implications for shifting species distributions (Hulme 2011), species abundance (Willis et al. 2008), local extinctions (Wien 2016), success of exotic species (Wolkovich et al. 2013), and intra-and interspecific interactions (Pardee et al. 2017, Robinson et al. 2017. Shifts in the herb community may enhance the light environment of early herb species, depending on the relative shift in the tree community that determines understory light availability (Heberling et al. 2019b).
The annual life cycle of herbs includes multiple "phenological phases," including emergence, leaf/shoot expansion, flowering/fruiting, senescence, and dormancy. These phases collectively determine survival and reproduction, and therefore the plant's evolutionary fitness. Each phase has a beginning date and an ending date, each a so-called "phenological event," as well as a duration. Important links and constraints may occur among events and phases so that it is critical to understand how they relate across the entire growing season (Wolkovich and Ettinger 2014). Post et al. (2008) argue that the focus should be on clusters of events, so-called "life history aggregates," to determine the extent of temporal coherence of events. Such an integrative approach has been applied to only a few herb species under experimental (Post et al. 2008, Jacques et al. 2015, Li et al. 2016 or observational conditions (Iler et al. 2019). Each phenological event of the annual life cycle is triggered by a unique combination of environmental factors interacting with intrinsic plant controls. Therefore, each event may not be equally responsive to climate change or may not change in the same direction, as found for migratory birds (Miles et al. 2017). Also, events are not independent (Donohue 2003, Post et al. 2008. Subsequent events can be predicted by earlier events, but also by duration between events (Ettinger et al. 2018). These so-called "legacy events" have been shown for a model plant species (Wilczek et al. 2010) and tree species (Fu et al. 2014, Keenan and Richardson 2015, Ettinger et al. 2018, Xie et al. 2018), but are unexplored in herb species. Because events (and phases) are integrated and sequential, opposing direction of change of a legacy event or amount of change in response to climate change may affect subsequent events. It is critical to study all events of the annual life cycle that represent different morphological/physiological functions for a plant, affecting its carbon gain (Kudo et al. 2008), nutrient uptake (Estiarte and Peñuelas 2015), and reproduction (Hedhly et al. 2009).
Multiple metrics characterize each phase. A given phase has a duration, conscribed by its initiation and ending dates. Some phases, e.g., flowering, also have a peak date. In studies of herbaceous flowering, only one phenological event, i.e., date of first flowering, is the common focus (Abu-Asab et al. 2001, Ellwood et al. 2013. Biases arise in analyses that rely on first-in-population event dates , Lind en 2011, and obscure the most biologically relevant metrics ). Rarely, date of peak flowering , Calinger et al. 2013 or duration of flowering Stenseth 1999, Bock et al. 2014) are studied. The direction or amount of change of first, peak and last date of a phenophase in relation to climate change may not be comparable (CaraDonna et al. 2014). Thus, a shift in first date may not indicate a change in duration. Changes in duration of phenophases in response to climate change may be as important as changes in phase initiation (Ettinger et al. 2018).
Three pieces of evidence would further the inference that a shift in phenology is associated with a shift in climate: (1) a change in event timing, (2) identification of a specific weather variable that predicts the event, and (3) a significant temporal trend in the event matches a significant temporal change in the weather variable. Even satisfying these three criteria are insufficient to reach a conclusion that the phenological shifts are caused by climate change. First, the difficulties of selecting the appropriate climate variable, e.g., measures of temperature, are well known (Kðrner and Hiltbrunner 2018, Chmura et al. 2019, Keenan et al. 2020. Temperature is a complex phenological cue. Additionally, although a climate factor can be identified as being a strong predictor of a phenological shift, such an approach ignores the possibility of multiple correlated factors, and especially causation at a mechanistic level. This observational study provides correlative evidence, and cannot claim causation. Experimental studies (Yoshie 2008, Foley et al. 2009), including molecular approaches (Wilczek et al. 2010), contribute to identifying causal mechanisms of phenology (Chmura et al. 2019).
Early studies focused on step 2 found that temperature, usually in the month preceding initiation of flowering, predicts the event (Sparks et al. 2000, Fitter andFitter 2002). A study with step 3 used a nonspecific weather variable, e.g., annual minimum temperature, to predict first flowering date (Abu-Asab et al. 2001). Not all components of monthly or seasonal temperatures have changed in the same pattern or to the same extent (Wuebbles et al. 2017). Most studies report some combination of two of the three steps (Bradley et al. 1999, Lambert et al. 2010, Beaubien and Hamann 2011. Some report the coincidence of the temporal pattern of phenological date and weather variable (Walther et al 2002, Amano et al. 2010), but do not indicate whether both temporal patterns changed significantly. A weather variable may predict a phenological event well, but either the phenological event and/or the weather variable may not be changing significantly through time. Menzel et al. (2006) in a study of European trees was the first to include all three steps, and Anderson et al. (2012) and Ellwood et al. (2013) provide explicitly such evidence for herb flowering.
This study characterizes the phenology of the entire annual cycle of 43 prominent herb species of a temperate deciduous forest from 1993 to 2017. Data for eight of these species from 1949 to 1964 collected in the same area by Charles Smith provide a longer perspective.
Overall, the study addresses the hypothesis that a change in date of a specific phenological event (beginning and ending) and duration of a phase in the annual life cycle is correlated with a specific component of climate change. To support the hypothesis requires demonstrating the three pieces of evidence developed above. The 25yr study based on observations by one observer has four special attributes. First, it is a community-level field study of herbaceous species of temperate deciduous forest. Second, it addresses temporal changes in all major phases during the entire annual life cycle. Third, it uses multiple metrics to address events (beginning and ending dates) and phases (duration). Fourth, it determines whether or not specific temporal changes of both temperature and precipitation, and earlier phenological events, match the temporal changes in the phenology of herbaceous species. By appraising how inter-annual changes in weather predict changes in all phenological components of the annual life cycle, this comprehensive long-term study provides the first integrative perspective of herbaceous species' associations with climate change.

Phenological data: Trelease woods 1993-2017
Observations of annual phenological events for 50 herbaceous species were made in Trelease Woods, a 24.5ha fragment of temperate deciduous forest, dominated by sugar maple near Urbana, Illinois, USA (40.084°N, 88.2404°W; see Augspurger and Buck 2017 for site details). Elevation varies by <5m. The 25 1-m 2 quadrats observed were ≥50 m from the edge in a 50 9 50 m grid in the north of Trelease Woods. The phenological status of each species in each quadrat was determined by the same observer (C. K. Augspurger) weekly in spring and summer, biweekly in fall, and biweekly to monthly in winter from 1993 to 2017 (Appendix S1: Text S1). For each species in each quadrat, dates of six phenological events were used to calculate five durations of phenophases (Table 1). Thus, measurements were recorded, and later analyzed, at the quadrat level. Quadrats differed in number and size of individuals; regardless, the first (or last) individual in the quadrat to begin or end a phase defined the event date for that species in that quadrat. No adjustments were made for this variation unless an interannual pattern in plant size was obvious, whereby the species was excluded from analysis.
Species varied from 1 to 25 quadrats and 6 to 25 yr of observation. Criteria for a species to be analyzed were ≥15 yr distributed relatively evenly over the 25 yr; a minimum or constancy in number of quadrats were not required. Added criteria were based on a species' natural history, including complex phenologies such as multiple leafing periods, and ability to apply consistent criteria of an event (Appendix S1: Text S1).
Species were analyzed together in seasonal groups (Appendix S2: Table S1). It was not desirable to group all species into a single analysis because, among species, the date for a given event ranges widely covering different months/seasons so that species are responding to different weather conditions. Also, it was not logistically practical to build climate models for 43 species separately for up to six events and five durations. The compromise solution was to categorize species into different months/seasons. For most events, species were grouped Notes: Change is represented as an advance or delay (for events) or shorter or longer (for phases). Also shown are the number of species analyzed and the percentage (in parentheses) of species advancing or shortening. Days of year (DOY) of the period covered by each group are provided. See Appendix S2: Table S1 for species in each seasonal group.
† No fixed effect was included in the top model. as early or late (Table 1, Appendix S2: Table S1). In addition, three categories were used for begin senescence and begin dormancy, or durations based upon these two events, only because of the longer period covered by these events over spring, summer, and autumn and there being natural divisions in the species' dates. Natural breaks with periods with no species active occurred for most events. Artificial breaks at the end of the month were used only for emergence and end expansion events. These artificial breaks were justified based on results from later analyses (see Weather variables: 1993Weather variables: -2017. Duration groups were based on the seasonality of the end event (final day) of the phase. The mean number of species in the seasonal groups in analyses of the six events and five durations was 10 (range 4-23 species). Eleven species were analyzed for all events and phases; an additional seven species were used in no analysis, although their data can be found in Fig. 1 and Dryad Data Repository (see Data Availability). For each combination of seasonal group and phenological metric, the number of records used in the analyses ranged from 2804 to 5483.
Phenological data : 1949-1964 Observations of the same six events of eight of the same species (genus in case of Viola spp.) were made at unknown frequency from 1949 to 1964 by Charles Smith, woods custodian. The location from 1949 to 1956 was noted as Brownfield Woods (5 yr) or Trelease Woods (6 yr) or both in the same year (3 yr); no location was given for 1957-1964. Prior to 19th century deforestation, the two woods were in the Big Grove. Today, the distance between the two fragments is %2 km. Descriptions by Smith are as follows: first appearance above ground (emergence by C. K. Augspurger), leaves full grown (end expansion), begin flowers (begin flowering), end flowers (end flowering), leaves begin to color (senescence), leaves all fallen/withered (dormancy Values represent interannual mean dates per phenological event; within a given year, the value for an event was calculated as the mean date across quadrats. Species are ordered by mean emergence date. Vegetative phenology is shown for leaf expansion (from emergence through full leaf expansion), mature leaf (from full leaf expansion until beginning of senescence), and senescence (from beginning of senescence until dormancy when no leaves remain). Flowering phenology is also shown (from beginning to end of flowering). Dotted light green lines indicate continuing output of new expanding leaves at least until flowering begins (six species) or leaf expansion for a second cohort of leaves (three species). Dotted dark green lines indicate the mature leaf period for the second cohort of leaves (four species). Dotted orange lines following dotted dark green lines indicate senescence of a second cohort of leaves (one species). Dotted orange lines following solid orange lines indicate species that retain leaves into January in some years (two species). Dotted orange lines following solid dark green lines indicate species with semi-senescent leaves that do not disappear until new growth in the subsequent year (five species). Aplectrum hyemale has leaves only from autumn to spring. Four species of winter annual species emerge as cotyledons in the autumn of the prior year and expand leaves in the subsequent year. Hydrophyllum appendiculatum is a biennial; H. appendiculatum yr 1 refers to its first year of leaves (with second cohort of leaves in some years) and yr 2 refers to its second year with new leaves, flowering, and senescence. Not all species were used in analysis of a given event (see Appendix S1: Text S1 and Appendix S2: Table S1 and Table S6); species names preceded by a dagger ( †) were used in no analysis. For 1949For -1964, one date per species per year was recorded for each event; no variation was noted. That date was assumed to represent mean values of a given event. For comparison, mean values for each event for each year for each of the same eight species for 1993-2017 were calculated. Differences between means for each of 44 species-events were determined by t tests.

Weather variables: 1993-2017
Daily temperature and precipitation data from 1993 to 2017 are from the Midwestern Regional Climate Center for the Champaign, Illinois Weather Station (3S), 8 km southwest of the study site (data available online). 2 The open area surrounding the station is similar to the edge of the study site, resulting in temperatures slightly above those of census plots in the forest interior.
Weather variables used in model selection to predict each event were based on literature for herb or woody plants (see Appendix S2: Table S2 for details). For the six events, temperature variables were (1) mean of the daily mean temperature ([(daily maximum temperature + daily minimum temperature)/2] for each day; then mean of month or season), (2) heating degree days (HDD; Σ (T b À T mean ), where T b is base temperature and T mean is daily mean temperature) at thresholds ranging from À4°C to 12°C, and (3) cooling degree days (CDD; equivalent to GDD, growing degree days; Σ (T mean À T b )) at thresholds from 0°to 12°C. Herb species emerge after over-wintering in the ground and may respond most directly to soil, not air, temperatures. Monthly mean soil temperatures at both 10.2 and 20.3 cm, available from only 2014 to the present, were highly correlated with monthly mean daily (air) temperature (r = 0.96, P < 0.0005). Precipitation variables were (1) total precipitation (P), (2) number of rainy days (RD), (3) number of heavy rain days with precipitation ≥20 mm (flood potential; HRD), and (4) number of weeks of drought (7 consecutive days with no precipitation; DR). Snow cover does not persist into February and March at the site. All phenology involves a sequence of phases and each phase may not be independent of the earlier phase(s) or its duration. Therefore, additional dates of emergence, end expansion, begin flowering, and/or begin senescence were added as variables ("legacy events") predicting some subsequent events. Durations from an earlier event were not used as explanatory variables. A specific combination of variables and months was used for each event, depending on its season (Appendix S1: Text S2, Appendix S2: Table S3).
Emergence and end expansion were analyzed first without seasonal groups. Including all species in the same analysis was problematic because their periods of emergence, March through April, and end expansion, April into June, cover a large range of dates when temperatures are changing greatly. Signs of the statistical artifact became evident in later steps (e.g., see Appendix S3: Fig. S1) and therefore all final analyses were done on seasonal groups.
Generally, the beginning of temperature-related variables preceded the event by 4-6 months and the precipitation-related variables preceded the event by 2 months; both temperature and precipitation variables continued through the end of the event for the seasonal group. For example, weather variables included in analyses for emergence and end expansion events stopped with the phenology of the species in those groups (e.g., April temperatures were used for emergence, late group, which included only species with the event date in April). For other events covering longer periods, months were included for which some species were no longer active (e.g., in the late group with date of first flowering covering 19 May-14 June, both May and June temperatures were used, even though some species had the first flowering date in May). In no analysis, however, was a later weather variable found to be the strongest predictor for events that occurred prior to that variable's period.

Statistical analysis
We used linear mixed-effects models to test for phenological changes through time and the weather variables that had the most predictive power of phenological events, as in Xie et al. (2018). Distinct approaches were used for the two sets of analyses. In all cases, analyses related to mixed-effects models were conducted using two packages in R version 3.5.1 (R Core Team 2018).
To address temporal trends in phenology, we tested 13 total phenological response variables in separate analyses for the dates of six events and 11 for the durations of five phases. Individual quadrat observations, in each year, were used in the analyses. For each phenological variable, we used model selection to find the top model at the seasonal-group-level with respect to fixed effects and random effects structure. All models included random effects terms for the species and quadrat, and the two were treated as crossed random effects (Baayen et al. 2008). Therefore, the analysis also provided species-specific temporal changes. We used a two-step approach in model selection that used separate information criteria for testing fixed and random effects (Xie et al. 2018): (1) AIC c (Akaike information criterion corrected for sample size) to determine the fixed effects structure, (2) conditional AIC, or cAIC (Greven and Kneib 2010), for the random effects structure. The cAIC is the preferred method for model selection based on random effects (Vaida and Blanchard 2005). When determining the fixed-effects structure, we compared random slope models with and without year as a fixed effect. We used the lme4 package (Bates et al. 2015) to construct mixed-effects models for this portion of the analysis. Models were initially constructed and compared using maximum likelihood. Once we arrived at the top model, final model parameters were determined using restricted maximum likelihood. Confidence intervals for the species-level coefficients of the top model were estimated with 10,000 simulations using the arm package (Gelman et al. 2016).
For analyses of the relationship between phenological events and weather variables (Appendix S1:Text S3), we tested the dates of six phenological events in 13 separate analyses based on seasonal groups of species (Table 4). Fixed effects consisted of weather variables, and in some cases phenological variables that preceded the response, the legacy effects. Weather variables described temperature and precipitation in individual months or seasons, and differed for each group of response variables (Appendix S2: Table S3). The random effect of species was included in all models. For the mixed-effects models describing the relationship between phenology and weather, we used the nlme package (Pinheiro et al. 2018) because it could account for heteroscedasticity in our model residuals.
We completed three main steps to select weather variables that were most predictive of an event.
(1) For each weather variable, individually determine the random and fixed effects structures for the top model predicting the phenological event. Eliminate from further consideration those weather variables that were not included as a fixed effect.
(2) Among correlated weather variables, choose the strongest predictor. Exclude weather variables not in the top model predicting the phenological event. Eliminate from further consideration those weather variables not included as a fixed effect, i.e., that are correlated with a stronger predictor. (3) Test different combinations of the remaining fixed effects, highlighting the top models and the weather variables they include. Further details of these three steps are in Appendix S1: Text S3.
We used two approaches to evaluate models; see Appendix S1: Text S4 for details.

Phenological patterns
Mean dates of each phenological event and each duration of phases differ greatly among herb species in this temperate deciduous forest community (Fig. 1). Although the categorization of herb species into "syndromes" is a useful heuristic device (Givnish 1987, Uemura 1994, some species of this study do not fit neatly in such discrete categories. Most durations of a given species cover multiple months and sometimes multiple seasons. Emergence by all but five species occurs between day of year (DOY) 65 and 120 in the spring. The duration of leaf expansion until leaf maturity can be as rapid as 9 d or prolonged to 211 d if the entire plant continues to grow from spring into summer/autumn. The period of mature leaves varies from very short for spring ephemeral species to year-round for "quasi-evergreen" species. Flowering generally is of shorter duration (15.8 AE 6.5 d; mean AE SD) and is most likely to occur in the spring. Senescence duration is shorter for spring species than for summer or autumn species.

Phenological change through time: 1993-2017
Dates of all phenological events changed from 1993 to 2017 for most seasonal groups (Table 1, Appendix S2:  Table S4). Both amount and direction of change in date varied at the seasonal-group-level among the six events (y-axis of Fig. 2 Fig. 2). For a given species, on average, 68.4% of events advanced and 74.4% of durations shortened; both are significantly more than expected at random (two-sided exact binomial test, P < 0.00001).
For a given event, the rate of change also varied among species (Fig. 2, Appendix S2: Table S6), sometimes showing divergent patterns in advance or delay (Table 1). In early species, for all events except begin dormancy, the degree to which a species shifted over time was correlated with that species' phenology, i.e., the top statistical model had random species slopes and intercepts that were correlated (Appendix S2: Table S4). For example, for early species, emergence showed less delay (y-axis in Fig. 2A, representing random species slopes) with later baseline emergence dates (x-axis in Fig. 2A, the random species intercepts). In late species, for all events except end flowering and begin dormancy, the degree to which a species shifted over time was related to its phenology, but patterns differed among events (e.g., Fig. 2B vs. Fig. 2F).
For early species, expansion and senescence durations were shortened; with greater shortening in species with longer initial durations (Fig. 3A, E, Appendix S2: Table S5). For late species, expansion, flowering, and senescence durations shortened (Fig. 3B, D, F). Emergence-senescence duration lengthened (Fig. 3H). Some phases showed no change in duration.
In model evaluations (see Appendix S1: Text S4), 10fold cross-validations of the 13 event-season combinations showed that model prediction uncertainties averaged 11.2 d (range 7.5-26.8, median 8.8), as measured by the root-mean-square error (RMSE). Correlations between predicted and observed dates yielded R 2 values that averaged 0.46 (range 0.15-0.65, median 0.5) (Appendix S3: Fig. S2). For the 11 phenophase-season combinations, 10-fold cross-validations showed that model prediction uncertainties averaged 12.7 d (range 4.4-27.8, median 11.4), and the R 2 of the correlation between observed dates and model predictions averaged 0.43 (range 0.24-0.64, median 0.43) (Appendix S3:   lays). Values on the x-axis indicate relative phenology (random species intercepts), where lesser values indicate phenologically earlier species. Solid vertical lines and dashed horizontal lines give the 95% bootstrap CI. Change in days represents the coefficient (random species slope, b) from linear mixed-effects models, while the model intercept is the random species intercept from those models, as in the equation date = b(year) + intercept, where the year 1993 is set to year = 0. Numbers refer to species identifiers (see Appendix S2: Table S1). The large red dot is the coefficient of the same variables from the seasonal group-level analysis (AE 2 SE; red vertical and horizontal lines) (Appendix S2: Table S5). For red dots with no vertical SE, the top model did not include year as a fixed effect, i.e., the event date did not change through time. Begin dormancy, early species was excluded because of a lack of random species slopes in the top model (Appendix S2: Table S4). Phenology in 1949-1964 vs. 1993-2017 Generally, events occurred earlier for 1993-2017 than 1949-1964 (29 of 44 species-events advanced; Table 2). A significant advance occurred in 23 cases, a significant delay in six cases, and no statistical difference (t test P > 0.05) in 15 cases. The initiating events of emergence and first flowering showed minor difference between the two periods (mean delay of 0.3 d/decade and 0.7 d/decade, respectively). The completion events of end expansion, end flowering, and begin dormancy were usually earlier in 1993-2017 (mean advance of 2.4-6.2 d/decade), resulting in shorter duration of expansion, flowering, and senescence. Laportea canadensis was the only species with consistently earlier phenology for 1949-1964 than 1993-2017. On average, the 44 species-events advanced a mean of 2.1 d/decade over the 48 yr between the median dates of the two data sets. This rate is similar to some slopes from the mixed models based on 1993-2017 phenologies (Table 1, Appendix 2: Table S5).

Weather factors predicting phenology: 1993-2017
For most phenological events, one weather variable had the strongest predictive power (Table 3) and was included consistently in the top statistical models (Table 3, Appendix S2: Table S7). The three exceptions were end flowering (late species) and begin dormancy (early and intermediate species). A temperature variable (mean daily mean from March-April or May, CDD for March, "early spring" or "late spring") was the strongest weather predictor in 10 of the 13 top models (Table 3). Higher spring temperatures predicted earlier emergence, end expansion, begin and end flowering, begin senescence for early species, and emergence, end expansion, begin flowering, and begin dormancy for late species. In contrast, higher September temperatures predicted later begin senescence (late species). Secondarily, greater HDD values predicted later emergence (late species) and begin/end flowering (early species).
Legacy of earlier events was the strongest predictor for three events and secondary for three other events (Table 3). Generally, later phenological event dates and earlier legacy events were positively correlated, except emergence was negatively correlated with begin senescence (late species). Precipitation variables played secondary/tertiary roles in predictive models (Table 3).
For each of the strongest weather predictors, random slopes fit to each species all had the same sign (i.e., the model suggested all species were correlated with the weather variable in the same direction). However, the amount of change varied among species (Fig. 4, Appendix S2: Table S8). In general, events advanced by 1.8-2.9 d/1°C increase, or 0.03-0.09 d/1-unit CCD increase (Table 4). Begin senescence of late species was delayed by 2.9 d/1°C increase (Table 4). In evaluating the 13 models with weather variables (see Appendix S1: Text S4), 10-fold cross-validations showed that model prediction uncertainties were about 8.2 d (range 4.2-15.5, media 6.6) and R 2 averaged 0.71 (range 0.56-0.84, median 0.72; Appendix S3: Fig. S4).

Random effect variance
For all model types, species intercept dominated random effects structure. In 13 models that tested for changes in event dates through time, and 11 models that tested for changes in phenophase durations, most variation was observed in species rather than quadrat intercepts. The SD of the quadrat intercept was at most 17.4% of the SD of the species intercept (mean 5.6%, median 2.7%), and the SD across 25 quadrats never exceeded 4.2 d.
In the models testing weather variables that predict phenological events, similar to Xie et al. (2018), random effects of species had much greater standard deviations than weather variables, in most cases with a greater than 100-fold difference (Appendix S2: Table S9).

Temporal match of phenology and weather: 1993-2017
The final step in support of the hypothesis that change in these phenological events is associated with climate change requires a demonstration of a significant temporal change in both event date and the weather predictor of that event during 1993-2017. Six of 13 event-season combinations met that requirement, viz. for early species, end expansion, begin flowering, end flowering, and begin dormancy; for late species, end expansion and begin senescence (Fig. 5, Appendix S2: Table S10). The significantly changing weather predictors were temperature variables associated with warmer "early spring" (a combination of March/April) or September.
The other seven events, emergence and begin senescence (early species), emergence, begin and end flowering (late species), and begin dormancy (intermediate species and late species), were not associated with climate change (Fig. 5, Appendix S2: Table S10). Five phenological variables did not change through time (Fig. 5B, F, I, L, M), although two were associated with a temperature variable that increased significantly (Fig. 5L, M, Appendix S2: Table S10), i.e., failure at step 1 from our three criteria. One event (end flowering, late species) did not have a strong predictive weather variable (Table 3), i.e., failure at step 2. Emergence (early species) advanced through time, and was predicted by a weather variable, but the weather variable did not change through time (Appendix S2: Table S10), i.e., failure at step 3. For the seven phenological events without support for a shift with climate change, only begin senescence (early species) had a secondary weather variable (March drought) that consistently appeared in the top models predicting the event (Appendix S2: Table S7).
However, this weather variable did not change significantly through time (Appendix S2: Table S10).

Temporal change in phenology
Temporal change in the phenology of this herbaceous community has been prevalent from 1993 to 2017. For early species, dates of four events are occurring earlier, while emergence has a slight delay and begin senescence is unchanged. For late species, end expansion and end flowering are advancing, while begin senescence is delayed. This high degree of phenotypic plasticity illustrates that the individual can match its activity to a changing environment (Anderson et al. 2012). However, change in some events may be maladaptive and a mechanism may have evolved to limit an event's plastic response to a changing environment. For example, emergence by March-emerging herb species delayed slightly over the 25-yr period. Some of those species are vulnerable to frost (Augspurger and Salk 2016). Emerging any earlier may be maladaptive because it would increase the frequency and duration of frost exposure (Augspurger 2013). The amount of change, advances of 3-4 d/decade for spring events of early species, are slightly higher than earlier estimates of 2.5 d/decade for leaf unfolding of European trees (Menzel et al. 2006). The late species' delay of 17.7 d/decade of begin senescence is markedly higher than the 1.3 d/decade delay in autumn leaf fall of European trees (Menzel et al. 2006). Flowering by herbs has advanced as much as 4.5 d/decade (Fitter and Fitter 2002), comparable to the 4 d/decade advance by this study's species, but all estimates are sensitive to which decade is compared (Amano et al. 2010).
Arising from these event changes, for early species, durations of three of five phases, expansion, senescence, and growing season (emergence-dormancy), became 2-4 d/decade shorter. For late species, spring phases were shorter (expansion duration by 8 d/decade and flowering duration by 2 d/decade), while autumn phases were longer (senescence duration by 9 d/decade and growing season (emergence-senescence) by 20 d/decade). Flowering duration of herb species on a British island shortened by 10 d/decade (Bock et al. 2014), but one-third of species bloomed longer by 13-19 d over 50 yr in Norway (Post and Stenseth 1999). No studies estimate changes in length of growing season of temperate herb communities, but European trees lengthened it by 10.8 d from 1959 to 1993 (Menzel and Fabian 1999).
This study indicates the value of using multiple metrics to describe each phase . Knowing only the beginning of a phase does not predict its duration. If both beginning and end changes are comparable, little to no change in duration occurs, as in our results for flowering duration (early species). Conversely, if beginning or end change in opposite directions, or if only one of the two changes, the duration of the phase  (Cleland et al. 2012), and the amount and rate of usage of stored reserves during dormancy (Risser and Cottam 1968).

Variation in change among events, phases, and species
The large variability among events and phases demonstrated by this study illustrates the importance of incorporating the entire annual life cycle in a phenological study. First, events and phases of this study's herb species did not change uniformly. For example, for both early and late species, advance was greater as the event occurred later in spring, i.e., end expansion and begin/end flowering had more change than emergence. For late species, expansion phase was shortened while senescence phase was lengthened. Second, within a given event, seasonality of a species affected the amount of change among species. For example, for early species, begin senescence did not change, while late species showed much delay. Third, within a given phase, duration of a phase affected the amount of change among species. For late species, those with longer flowering duration and growing season durations had a greater change in duration. These patterns related to seasonality of event and duration of phases are related to weather patterns (see Weather predicts phenology). Only short-term studies, using either experimental warming (Price and Waser 1998, Cleland et al. 2006, Post et al. 2008 or observations (Li et al. 2016, Heberling et al. 2019a, have been done on multiple phenological phases of multiple herb species. Although past long-term observational studies found greater advancement of species with early than late spring flowering (e.g., Fitter andFitter 2002, Miller-Rushing and, this pattern cannot be assumed to be predictive of other events in the annual life cycle. The lack of a consistent pattern of change among events and durations limits any simple, comprehensive summary of phenological change of herb species. All events were not expected to advance. Spring warmth advancing flowering has been the primary focus of previous studies, but senescence may delay with warming. In fact, species did not show a temporal coherence of shifts within the annual life cycle. They had limited aggregate life history responses (Post et al. 2008). On average, a species had 68.4% of events or phases change. Some species had five of six events advance or all phases shorten, while other species had only two in six events advance or two in five phases shorten (Appendix S2: Table S6). Only one aggregate life history response was obvious. Begin and end flowering advanced for all species. The legacy event of begin flowering was a strong predictor of end flowering in the climate models. Additionally, the two growing season length metrics, summarized as emergence-senescence or emergence-dormancy, responded in the same direction. Senescence was a strong predictor of dormancy in the climate models. Except for begin senescence of late species, no legacy event from spring extended its influence into autumn. In general, the predictive legacy event was in the immediately preceding phase and in the same season. Apparently, the developmental sequence initiated by early weather was uninterrupted by a secondary, contrasting weather driver or the two phases have an endogenous dependency (Lechowicz 1995). No comparative studies of legacy events are available for herbs. For temperate trees, dates of senescence were related to dates of leaf out in one study (Keenan and Richardson 2015), but not in another (Panchen et al. 2015). Ettinger et al. (2018) found that early events of trees weakly constrain most later events. Generally, while patterns relating to Notes: Analyses were completed for each of six phenological events, with different seasonal groups (13 total models). See Appendix S2: Table S2 and Table S3 for descriptions of each weather variable. Also shown is the number of additional occurrences of each variable in the next four top models. See Appendix S2: Table S7 for the next four top models.  (Table 3) and random legacy events are beginning to be found, the underlying mechanism is yet to be explored.
All events may not be equally plastic to environmental change. Only four of six events showed a large advance in this study. In an earlier experimental study of six of this study's spring ephemeral species, six cohorts were forced experimentally to emerge over a period from January through June, but leaf senescence occurred in June for all cohorts (Augspurger and Salk 2016). Events of the entire annual life cycle were not able to adjust to a change in one earlier event. Any lack of aggregate, coordinated responses throughout the annual life cycle has consequences because the individual has to adjust to a combination of advances and delays of event dates, shortening or lengthening of phases, and new timing of subsequent events and phases. Any shift may result in a specific phase not matching favorable conditions for survival, growth, and reproduction. These complex changes and subsequent adjustments may not be optimal for fitness (Post et al. 2008).
The large variability among species demonstrated by this study illustrates the problems inherent in restricting analyses to the community-or seasonal-group level.
Finding variation among herb species' phenological change was not surprising, given their wide range of syndromes, principally related to length and timing of growing season, i.e., when emergence and senescence (or dormancy) occur. Seasonality of the event of a given species impacted greatly the change in direction and magnitude of change. While for end expansion and begin/end of flowering, almost all species showed an advance, emergence and begin senescence had delays as well as advances among species. For future exploration of effects of evolutionary relationships on phenological changes of species, the addition of "family" to the hierarchical structure of the fixed effect model could be used. This interspecific variation among herb species in the amount and direction of change in phenology has consequences for the species composition and relative abundance of species in a community (Roberts et al. 2015).

Weather predicts phenology
Among the limited variables tested, a temperature variable was the strongest predictor of phenology for slopes varied among species (Appendix S2: Table S8; for events with slopes that did not vary among species, slopes are in Table 4). Change in date per weather variable unit is the random-species slope for the top model describing the relationship (Table 3,  Appendix S2: Table S9), while relative phenology date is the random-species intercept from the model. Smaller model intercept (xaxis) values correspond with earlier phenological events. Numbers refer to species identifiers (see Appendix S2: Table S1 for species). The large red dot is the coefficient of the same variables from the seasonal group-level analysis (AE2 SE; red vertical and horizontal lines). CDD, cooling degree days; subscript numbers indicate the base temperature (T b ) in°C. (FIG. 4. Continued) TABLE 4. For each phenological event, a summary of the seasonal group-level shift (number of days) per unit change in the top weather variable (see Table 3). Notes: Also shown are the number (percentage in parentheses) of species with a negative slope in the weather model, indicating an advance with an increase in the weather variable. See Appendix S2: Table S1 for species in each seasonal group. Number of species in each group and percentage of species (in parentheses) with a negative slope.

Event and group
† Date of begin flowering, and not any weather variable, performed well as phenological explanatory variable. ‡ The strongest weather predictor, but date of begin senescence performed stronger as phenological explanatory variable.  Table S7), the variable chosen had the highest absolute t value most events. Warmer spring temperatures predicted earlier phenological events. A temperature variable also was the second-strongest predictor for the three exceptions that had prior phenological events, so-called "legacy events," as the strongest predictor. For spring events, the mean daily mean temperature of the preceding month (May or March) and/or March/April CDD (of different thresholds) predicted the event. Those temperature variables are highly correlated with other temperature variables, including mean daily temperatures and maximum and minimum temperatures, as well as soil temperatures. Accumulating warmth from rising spring temperatures is well known to force plant development, including flowering of temperate herbs (Fitter andFitter 2002, Miller-Rushing and and bud burst and leaf expansion of temperate deciduous trees (Polgar and Primack 2011). For begin/end flowering of spring species, HDD was a secondary predictor, perhaps indicating an additional chilling requirement prior to spring warmth for these events (Michaels and Amasino 2008).
For autumn events, higher temperature of the preceding month predicted a delay in date of begin senescence. Spring CDD was a surprising secondary predictor of begin dormancy (late species). Senescence date was the strongest predictor, as for trees (Xie et al. 2018). This result may indicate that the analysis did not uncover the actual weather predictor or there is no weather predictor. Low samples of intermediate and late species for begin senescence and begin dormancy may have limited detection of their response to weather. Trees show a delay of autumn events although the pattern is less strong than the advance of leaf expansion in spring (Menzel and Fabian 1999). Both temperature and precipitation affect leaf senescence and drop of trees (Archetti et al. 2013). Warm temperatures and drought delay autumn phenology of trees, while heavy rainfall and heat stress advance leaf senescence and drop (Xie et al. 2018). Precipitation variables of spring only, including monthly precipitation, number of rain days, number of flood days, and weeks of drought, were of secondary or tertiary importance as predictors of herb phenology.

Thermal sensitivities of events and species
Any change in a species' phenology depends on both the amount of temperature change and the species' thermal response (Shen et al. 2014). Although thermal sensitivities were calculated from the study's results, any relationship found between temperature and phenological change is purely correlational. Also, the weather predictor used to determine thermal sensitivities, e.g., integrated CDD in this study, has associated uncertainty and changes in variance depending on length of integration (Keenan et al. 2020). Thus it is only a proxy to the real integrated temperature to which plants respond.
Advances of the events averaged 2.3 d per°C increase or 0.06 d per CDD increase (Table 4). These estimates are seasonal-group wide, while individual species differed greatly in these values (Appendix S2: Table S8). These mean thermal sensitivities are comparable to previous studies of herb flowering of 2.4 d advance per°C (Calinger et al. 2013) and spring/summer events by European trees of 2.5 d advance per°C (Menzel et al. 2006). Late species begin senescence delayed by 2.9 d per°C, more than the 1.0 d per°C delay of leaf coloring/fall by European trees (Menzel et al. 2006). Events differed in their sensitivities. Begin senescence and begin dormancy of early species were less sensitive than other events.
Thermal sensitivities varied two-fold among species, e.g., for emergence, from 1.2 to 2.6 d per 1°C and for end expansion, from 0.048 to 0.108 d per CDD. Flowering of herbs in general shows a sensitivity of 5 d per 1°C increase (Amano et al. 2010), but can range among species from 2 to 10 d per 1°C (Sparks et al. 2000) or 3 to 6 d per 1°C . This study's begin and end flowering of early species each advanced 0.06 d per CDD (the equivalent of about 1.9 d per 1°C increase). Therefore, species shift the same event at different rates and a given species shifts its events at different rates, explaining some of the complexities of temporal change throughout the annual life cycle. Mechanisms underlying phenology are likely different from event to event and species to species. Studies of differences in flowering responses find that some species may be controlled by factors other than temperature or are no longer experiencing the appropriate temperature pattern to cause vernalization to occur (Cook et al. 2012). Some species may simply have a very low sensitivity to temperature differences or lack adaptive plasticity to respond to such change.

Association of phenological events with climate change
This community-wide, 25-yr study of the annual phenology of temperate herbaceous species provides some strong support for the hypothesis that phenological change is correlated with climate change. The three required interconnected lines of evidence were met for (FIG. 5. Continued) (Table 3). An inverted scale for the event date (left y-axis) was used for negatively correlated variables (all panels except J). For weather variables, regression lines represent the least-squares regression through time; solid lines (P < 0.05) and dashed lines nonsignificant results (Appendix S2: Table S10). For phenological events, regression lines were created from the parameters in the change-through-time mixed-effects model (Appendix S2: Table S5); solid lines represent cases in which year was included in the top mixed-effects model; dashed lines represent cases in which year was not included in the top model (Table 1; Appendix S2: Table S5). An asterisk (*) indicates that both phenology and weather variables have significant trends through time (P ≤ 0.05), indicating support for the hypothesis that phenology has shifted with climate change. All dates are day of year. some event-season combinations, but not all. First, event dates change substantially for end expansion (early and late species), begin and end flowering (early species), begin senescence (late species), and begin dormancy (early species) (Table 1; Fig. 2). Second, a specific weather variable predicts well the date of each of these event-season combinations (Table 3). Third, and importantly, the temporal pattern of phenological change tracks the temporal change in the specific weather predictor for each of event-season combinations ( Fig. 5; Appendix S2: Table S10). This support is bolstered by phenological advances by seven of eight of the same species in the same area in 1993-2017 relative to 1949-1964 (Table 2). Therefore, at least this subset of temperate herb species has been responding to climate change for at least 50 yr. Comparable long-term studies of multiple events of herb species are not available for comparison.

Lack of concordance in climate and phenology
Seven of 13 event-season combinations do not show evidence of a shift in phenology that is associated with climate change (Fig. 5; Appendix S2: Table S10). All of those combinations, except end flowering in late species, have a strong weather predictor. However, the lack of concordance arises because either the date of the event and/or the weather variable does not change with a significant trend through 1993-2017. Temperature changes over the past decades in the contiguous United States (Liu 2017, Wuebbles et al. 2017) and in Illinois during the 20th century (Robeson 2002) have not been the same for each month of the year. A nonsignificant temporal change in phenology in our seasonal group-level analysis likely obscures underlying matches for individual species. Species collectively generating the seasonal-group-level response likely differ in their temperature requirements to initiate a given event, in their thermal sensitivity, and/ or their specific mechanisms that drive their phenological responses.

Limitations
Even meeting all three required pieces of evidence allows a claim of only an association and not a causal relationship between phenology and climate change. The results simply show a correlation between a phenology variable and a weather variable (and perhaps the wrong weather variable). Alternative hypotheses were not tested, e.g., deer herbivory and invasion by shrubs affecting plant size and abundance (Christopher et al. 2014), indirectly affecting phenology, or indirect effects on phenology from species interactions with mycorrhizal symbionts (Allsup andLankau 2019, Bennett andClassen 2020). Environmental vs. organismal explanations of phenological shifts need to be distinguished (Chmura et al. 2019). The predictor variables used in this (and most other) studies do not match any theory of the underlying drivers of plant phenology (Kðrner and Hiltbrunner 2018).
Integrating temperatures over a period of time to generate a predictor variable does not guarantee that a plant is doing the same integration of temperatures and over the same period (Keenan et al. 2020). Also, reducing analyses to simple predictors to cover multiple species included some species whose phenology does not match exactly the summation period. Restricting the analysis to only the top three predictor variables may have limited the ability to identify the correct weather variable(s). It is likely that multiple variables interact to determine phenology for many species. For example, two variables, perhaps chilling, then warmth (Michaels andAmasino 2008, Cook et al. 2012) or photoperiod and warmth (Bassler and Kðrner 2012, Flynn and, may interact to determine emergence dates. However, we did not test interactions. Finally, the use of linear regression to test for trends may limit a true test for an association between weather and phenology. Relationships may be nonlinear, and thus poorly suited to our linear-mixedeffects model approach (Jochner et al. 2016(Jochner et al. , 2018. Weather has large variability from year to year and climate change is not strictly linear, given the amount of "noise" generated by weather variability (Keenan et al. 2020).
Accounting for the variation in responses among events, phases, and species will require gaining a better understanding of the underlying mechanisms controlling phenology ). More studies for individual species are needed, e.g., the genetic and physiological underpinnings of the plant's response to its environment, its thermal sensitivity, and its integration of various events within the annual life cycle (Wilczek et al. 2010). Common gardens with experimental manipulation provide direct tests of possible mechanisms underlying phenological change (Berend et al. 2019). Such species-specific studies of herb species are slowly accumulating (Yoshie and Fukuda 1994, Lapointe 2001, Rothstein and Zak 2001, Jochum et al. 2007, Yoshie 2008, Estiarte and Peñuelas 2015. Any community-level study of many species necessitates compromises of balancing sample sizes of year, quadrats, and density of individuals. Ideally, a long-term study includes many representative species, as in this study. This forest's loss of individual herbs over the 25-yr period (Augspurger and Buck 2017) restricted the availability of sufficient observations to include some species in the analysis. Additionally, shifting population abundance and plant size in our study may have contributed to some perceived changes. Smaller populations can show a delay in the first date of an event . Smaller non-reproductive individuals can have earlier senescence than larger reproductive individuals (C. Augspurger, personal observation).

Future directions
At first glance, this study suggests that an advance in phenology in spring has a positive effect on the light environment of these understory species, some that expand leaves before canopy closure. However, a given herb species' gain by an advance in phenology can be evaluated only by knowing the concurrent shifts of canopy trees that regulate understory light levels (Heberling et al. 2019b). Carbon gain in spring before canopy leaf-out varies greatly among spring-emerging species and can be considerable in autumn after canopy leaf fall (Heberling et al. 2019a). A concurrent study of phenology of the canopy tree species at our study site allows an appraisal of the effect of change in tree phenology on each herb species' light environment (C. K. Augspurger and C. F. Salk, unpublished manuscript).

CONCLUSIONS
This first community-wide phenological study of annual life cycle events of temperate herb species showed much temporal change over the past 25 yr, and for some species for 50 yr, but with great variation among events, phases, and species. Simple generalities did not arise. Study of one event is clearly not predictive of the entire annual cycle. Documentation of all events and phases, and multiple metrics of each event, are necessary to describe temporal change. Not all events are correlated with climate change to the same extent. The study outlines the complexity of responses occurring at both community and species levels. It underscores the need for indepth studies of genetic and physiological underpinnings of phenology and thermal sensitivity of individual species, as well as comparative community-wide studies of more forests and biomes.