Royal Society Publishing

Forecasting phenology under global warming

Inés Ibáñez, Richard B. Primack, Abraham J. Miller-Rushing, Elizabeth Ellwood, Hiroyoshi Higuchi, Sang Don Lee, Hiromi Kobori, John A. Silander

Abstract

As a consequence of warming temperatures around the world, spring and autumn phenologies have been shifting, with corresponding changes in the length of the growing season. Our understanding of the spatial and interspecific variation of these changes, however, is limited. Not all species are responding similarly, and there is significant spatial variation in responses even within species. This spatial and interspecific variation complicates efforts to predict phenological responses to ongoing climate change, but must be incorporated in order to build reliable forecasts. Here, we use a long-term dataset (1953–2005) of plant phenological events in spring (flowering and leaf out) and autumn (leaf colouring and leaf fall) throughout Japan and South Korea to build forecasts that account for these sources of variability. Specifically, we used hierarchical models to incorporate the spatial variability in phenological responses to temperature to then forecast species' overall and site-specific responses to global warming. We found that for most species, spring phenology is advancing and autumn phenology is getting later, with the timing of events changing more quickly in autumn compared with the spring. Temporal trends and phenological responses to temperature in East Asia contrasted with results from comparable studies in Europe, where spring events are changing more rapidly than are autumn events. Our results emphasize the need to study multiple species at many sites to understand and forecast regional changes in phenology.

1. Introduction

In temperate ecosystems throughout the world, the timing of phenological events is shifting, and these shifts have been linked to recent global warming (Parmesan & Yohe 2003; Root et al. 2003; Menzel et al. 2006). Changes in phenology could have substantial repercussions for conservation of natural systems, potentially creating ecological mismatches between interacting species (Both et al. 2006; Post & Forchhammer 2008; Thomson 2010) or between species and their abiotic environment (Inouye 2008). Such mismatches could in turn alter the dynamics, functioning and composition of many ecosystems (Winder & Schindler 2004; Willis et al. 2008; Miller-Rushing et al. 2010). While some species are likely to be harmed by changes in climate and phenology, others may remain unaffected or be favoured under changing conditions. Discerning the likelihood of these different outcomes will depend on reliable forecasts of species' phenologies.

One of the major effects of shifts in plant phenology will be their impact on the length of the growing season. The length of the growing season has great ecological and biogeochemical significance for biological communities. It controls the number of days that food is available for a wide range of animals, especially insects. It is also an important factor in determining how much water is returned to the atmosphere via evapotranspiration, how much carbon is sequestered in new growth and how much nitrogen and other nutrients is absorbed from the soil (Norby et al. 2003; Delpierre et al. 2009; Richardson et al. 2010). Documented consequences of changes in growing season length are only just emerging, largely because of a general lack of data on growing season length and other key ecological variables (Menzel & Fabian 1999; Schwartz et al. 2006). At the biome level, extended growing seasons owing to global warming have been connected with increased gross primary productivity in forests (Richardson et al. 2009) and zooplankton productivity in lakes (Shuter & Ing 1997). Extended growing seasons have also lead to increases in the number of broods that some migratory birds have in a given year (Monroe et al. 2009), and will probably increase the number of generations many insects, including some important forest pests, go through in a season (Tobin et al. 2008; Jönsson et al. 2009).

Although records of phenological events have thus far been used to document species responses to recent global warming (e.g. Dunn & Winkler 1999; Fitter & Fitter 2002; Menzel et al. 2006), they can also be used to forecast phenological trends under future climate scenarios and at locations for which we do not have any records (Primack et al. 2009). However, such forecasting is difficult for at least two reasons. First, recent analyses have shown that species' responses to temperature can vary significantly among sites, suggesting that phenological changes at one location may not always be good indicators of changes at other locations (Sparks et al. 2005; Menzel et al. 2006; Gordo & Sanz 2009; Primack et al. 2009). This spatial variability poses particular challenges to forecasting at locations for which few or no data are available. Second, even when records exist for locations of interest, forecasts will typically need to extrapolate beyond the range of temperatures observed at a particular site.

When accounted for in models, these dual challenges of spatial variability and limited climatic range of observations will introduce greater uncertainty to forecasts, particularly for locations with few data and large predicted climatic changes. Process-based phenological models (e.g. Chuine 2000) have aided efforts to forecast phenology by incorporating a mechanistic description of species responses to temperature variability. Nevertheless, they do not directly address these two issues.

Here, we implement an analytical framework, a hierarchical approach, that can help to overcome those difficulties and greatly improve forecasts of phenological changes. We first calculate species-specific parameters that allow us to make realistic inferences across the range of each species. Second, we use those data, which encompass a wide range of climate conditions, to forecast the species' phenological responses to future climate change. Our forecasting potential will be somewhat limited because we are unable to predict species' responses to novel climates (Williams & Jackson 2007), and because predictions into areas with no data will be broad, reflecting the range of responses found in the whole dataset. Yet, forecasts following the hierarchical method will be realistic, because they will reflect the natural variability inherent in the data.

2. Material and Methods

(a) Phenology and climate data

Since 1953, the Japanese Meteorological Agency has been gathering data on over 120 phenological events (the timing of life-cycle events) for both plants and animals in the grounds of 102 of their weather stations (data available from Japan Meteorological Agency, http://www.jmbsc.or.jp/english/index-e.html). These weather stations are located across the latitudinal range of Japan, from northern Hokkaido (latitude 45°24.9′ N) to the southern islands (latitude 24°20.2′ N), with a corresponding gradient in annual temperatures from subtropical to boreal (figure 1). Within that latitudinal gradient, the Weather Service of South Korea has been gathering similar data on 20 phenological events at 74 weather stations across the country, with some observations dating from 1921 (data available from Korea Meteorological Administration, http://www.web.kma.go.kr/edu/unv/agricultural/seasonob/1173374_1389.html; figure 1).

Figure 1.

Locations (dots) of the meteorological stations in Japan and South Korea where phenological records have been collected.

Climate analyses have shown that this region has warmed in recent years (Primack et al. 2009), and models forecast that local temperatures will continue to rise in the future owing to both global warming and increasing urbanization in the region (IPCC 2007). In particular, winter temperatures warmed by an average of 1.2°C from 1953 to 2005, with higher rates of warming in areas with particularly high human populations (Primack et al. 2009). There was, however, no latitudinal trend in warming—higher latitudes did not warm more than lower latitudes (Primack et al. 2009).

The phenological observations were made by agency employees according to carefully defined sampling protocols that have remained constant for the duration of the study period (JMA 1985). The protocols directed that observers note the phenological activity of one individual plant for each species close to each weather station. From this dataset, we analysed records of plant phenological events in spring (nine species for 11 phenological events) and autumn (four species for five events; table 1). For spring, we examined first flowering dates for Prunus yedoensis (Yoshino cherry, a tree), Forsythia koreana (forsythia, a shrub), Prunus persica (peach, a tree), Rhododendron kaempferi (Kaempferi azalea, a shrub), Rhododendron mucronulatum (Korean rhododendron, a shrub), Viola mandshurica (Manchurian violet, a forb) and Zoysia japonica (Japanese lawn grass), and leaf-out dates for F. koreana, Ginkgo biloba (maidenhair tree), Morus bombycis (mulberry, a tree) and R. mucronulatum. Autumn phenology data were recorded for change in leaf colour for Acer palmatum (Japanese maple, a tree) and G. biloba, and for leaf fall in G. biloba, M. bombycis and P. yedoensis. All species were observed in Japan, except F. koreana and R. mucronulatum for which we have data only from South Korea. Additionally, P. yedoensis first flowering and G. biloba leaf fall were recorded in both Japan and South Korea. For three species, G. biloba, M. bombycis and P. yedoensis, we have both spring and autumn phenological observations at several sites in Japan. For those species, we examined the length of their growing seasons, defined as the number of days from leaf out to leaf fall for G. biloba and M. bombycis and as the number of days from first flowering to leaf fall for P. yedoensis. Leaf-out data for Japan were not available for this species, but in the South Korean records, first flowering date and leaf out are strongly correlated (ρ: 0.75). Because we were not concerned about the actual dates but rather the magnitude and direction of the response to changing temperatures, we considered first flowering date to be a good substitute to indicate changes in the beginning of the growing season.

View this table:
Table 1.

Time-series trends, number of positive and negative slope coefficient values for the regressions of phenology versus year at each location. The number of sites with slope value statistically different from zero is indicated in parentheses. Events: LOD, leafing-out date; FFT, first flowering time; LCD, leaf colouring date; LFA, leaf fall. N, total number of observations; S, number of sites; T, period of time during the year for which mean monthly temperatures were used in the analyses.

Many other investigations have previously used these data to document how spring and autumn phenological events and growing season length have changed over time (Kai et al. 1996; Chen 2003; Matsumoto et al. 2003; Ho et al. 2006; Doi 2007, 2008; Doi & Katano 2008; Doi & Takahashi 2008; Doi et al. 2008). However, none of these studies has attempted to use these results to forecast future changes in phenology and their likely impacts as we do here.

Given the magnitude of the dataset (11–150 sites and over five decades of sampling) and the standardized sampling protocol, we are confident our analysis provided us with a reliable estimate of each species' response to interannual variation in temperature and other factors.

(b) A model for plant phenology

We began our analysis by examining trends in the phenologies of the studied species at each location for which data were available. For these exploratory analyses, we ran simple linear regressions (phenology versus year) at each location using R's ‘lm’ function (R Development Core Team 2008).

We next explored the periods of the year that best explained the phenology of each species along their distributional ranges (correlations not shown; periods used in the analysis (highest correlations) are specified in table 1). We used monthly mean temperatures for these analyses, as these were the best temperature data available to us (i.e. we did not have access to daily data). Previous studies have shown that temperature accounts for a substantial portion of the interannual variation (more than other climate factors such as precipitation) in plant phenology in temperate regions (Menzel et al. 2006; Prieto et al. 2009), and in Japan specifically (Doi 2007). For the final analysis, we used temperatures from the period that best predicted the entire dataset for each species, even if the best period differed from site to site. This approach provided one consistent parameter across all sites and years that allowed us to estimate the overall response of the species to the temperature gradient. Parameters defining such responses could then be used to forecast species phenology at locations for which we do not have any data.

We analysed the phenology data as a function of temperature with several alternative models (table 2). Here, we describe the model that worked best for most species, model 6. In this model, phenology at site s in year t is estimated as a linear function of temperature. Year random effects, time, are added to the Gaussian error (ɛ) in a hierarchical (or multilevel) model where intercept (α) and slope (β) coefficients are estimated for each location from common prior distributions (see below):Embedded Image

View this table:
Table 2.

Models evaluated in the analysis of phenological responses to temperature. Linear, exponential and logarithmic models were chosen for comparisons as these were the best describing the data. a, α, intercepts; b, β, coefficients for temperature; ɛ, error term; γ, coefficient for latitude; time, year random effects; location, site random effects; w, spatially explicit random effects; s, subindex for site; t, subindex for year.

This hierarchical approach, in which we fitted a different line for each site (parameters αs and βs), allowed us to combine data from all the sites to estimate the species' overall response to temperature, while allowing for variation in this response among sites. Parameter estimation for one site was informed by data from other sites through the hyperparameter link (figure 2). This hierarchical structure was critical for making inferences at locations for which we do not have data and under future climate scenarios. It also allowed us to work with different datasets, accounting for the differences between them (e.g. countries and data collectors; Cressie et al. 2009; Ibáñez et al. 2009).

Figure 2.

Hierarchical structure of the analysis. Site-specific parameters are informed by data from other sites through the hyperparameters. Hyperparameters are estimated at the species level, while site-specific parameters are estimated for each species at each site.

To estimate the large number of parameters involved in this hierarchical model, we used Bayesian methods, which are well suited for hierarchical models (Gelman & Hill 2007). We estimated αs and βs from prior distributions αs ∼ Normal(a, Embedded Image) and βs ∼ Normal(b, Embedded Image), respectively. The prior parameters of these distributions, hyperparameters (a, Embedded Image, b, Embedded Image), were then estimated from distributions with uninformative parameter values: a ∼ Normal(0, 1000), σa ∼ Uniform(0, 1000), b ∼ Normal(0, 1000) and σb ∼ Uniform(0, 1000). The random effects associated with each year, timet, were estimated in a similar way: timet Normal(timemean, Embedded Image), timemean Normal(0, 10 000) and σtime ∼ Uniform(0, 1000). The error term, ɛst, follows a normal distribution, ɛst ∼ Normal(0, Embedded Image) and 1/Embedded Image ∼ Gamma(0.003, 1). These prior parameter values of the Gamma distribution constrain the variability of the phenological response to fall within the range of the observed data. The time-associated random effects took into account year-to-year variability that was not reflected in the temperature data. Such variability could be caused by other climate factors such as precipitation, wind or sampling factors. The individual random effects, or error term, accounted for any unexplained variability at each site and year, e.g. changes in observed individuals or in site conditions other than temperature. Parameters, a, Embedded Image, b, Embedded Image, timemean, Embedded Image and Embedded Image, estimated at the species level, were then used to compare species' overall responses to increasing temperature (figures 3 and 4, black lines).

Figure 3.

Spring phenology records (circles) and predicted overall species spring phenology response to temperature (black lines, mean (parameters a and b) and 95% predictive interval). Colours represent data (dots) and predicted local responses (lines (parameters αs and βs)) at two locations.

The Bayesian approach was also used to accommodate the uncertainty in the temperature data, because the monthly records available to use were a coarse approximation of the actual conditions experienced by individuals. By including temperature in the model as a latent variable that needed to be estimated from observed temperatures, the uncertainty in temperature experienced by individuals was included in the model (Clark et al. 2003). Temperature values used in the analysis were estimated as a function of the monthly records for the selected period (tmr, temperature monthly records) within the range of variability (Embedded Image: temperature variance) observed in the data; then, temperaturest Normal(tmrst, Embedded Image).

Posterior densities of the parameters were obtained by Gibbs sampling (Geman & Geman 1984) using OpenBUGS 1.4 (Thomas et al. 2006). Simulations were run for 75 000 iterations. Convergence was assessed from multiple chains with different initial conditions and required from 1000 to 10 000 iterations. Pre-convergence ‘burn-in’ iterations were discarded. Model selection was based on deviance information criterion (DIC; Spiegelhalter et al. 2006; table 3). The effective number of parameters was approximated by subtracting the deviance of the posterior means of the parameters from the posterior mean of the deviance. Adding this value to the posterior mean deviance gives a DIC for comparing models, where the best predictor of the data is the model with the lowest DIC. Overall predicted phenology as a function of temperature (mean and 95% credible interval) was estimated for each species (figures 3 and 4, black lines), together with the predicted response at two different locations (figures 3 and 4, red and blue lines). Plots of predicted versus observed phenology were also used to evaluate the fit of the model.

View this table:
Table 3.

Species analysed for spring and autumn phenological responses to temperature. Events: LOD, leafing-out date; FFT, first flowering time; LCD, leaf colouring date; LFA, leaf fall. Model columns show DIC values, bold numbers indicate the best model fitting the data (lower DIC).

3. Results

(a) Temporal trends

Regression analyses of phenology date versus year show the disparities among species and among locations within a species (table 1). For the period for which we have data, 1953–2005, seven of the 11 spring events analysed showed trends towards earlier spring phenology at most locations. Overall, autumn phenologies are occurring later over time. In general, delays in autumn phenology are more common than are advances in spring phenology (table 1). Examining records from the species for which we have both spring and autumn data, we find that P. yedoensis spring flowering is occurring earlier at 92 per cent of the 150 locations for which we have data, whereas autumn leaf fall is occurring later at all 11 locations. Morus bombycis spring leaf out is occurring earlier at 64 per cent of sites compared with 77 per cent of sites where autumn leaf fall is occurring later. Gingko biloba spring leaf out is occurring earlier at 82 per cent of locations, but autumn leaf fall and leaf yellowing are occurring later at 85 and 79 per cent of locations, respectively. The results are similar if we confine the comparison to locations with statistically significant trends (table 1). For these same three species, the autumn events are changing more quickly than are the spring events—i.e. the slopes of autumn events tend to have greater absolute values than do the spring events, as determined by a Wilcoxon's signed rank test (p < 0.001) (for slope values, see electronic supplementary material, table SA).

The combination of earlier springs and later autumns results in a longer growing season for all three species for which we have data on both spring and autumn events, although there is variability from site to site (table 1). There are no significant correlations between latitude and the slope associated with spring phenology, autumn phenology and the length of the growing season for P. yedoensis and M. bombycis (for correlation values, see electronic supplementary material, table SB). However, for G. biloba, there is a significant positive correlation between latitude and spring phenology and a significant negative correlation between latitude and autumn phenology, indicating that growing seasons are lengthening more rapidly at higher latitudes (electronic supplementary material, table SB).

(b) Model for plant phenology

We examined several models (table 2) to analyse changes in spring and autumn phenology as a response to changes in temperature. Model 6, a mixed model with time (year) random effects that estimated phenology as a linear function of temperature, fitted the data best for eight of the 11 spring events analysed and for four of the five autumn datasets (table 3). To facilitate comparisons among species, we only report results of model 6.

(i) Spring phenology

Results from the analysis of the spring phenological data show an overwhelming response to temperature, with all events tending to occur earlier at warmer temperatures (table 4, parameter b; figure 3). Ten out of 11 coefficients were significantly different from zero. The range of responses varied from 7.7 (s.d. = 1.17) days earlier per degree Celsius phenology for P. persica to a non-significant value of 0.35 (0.52) days earlier per degree Celsius phenology for R. kaempferi. Species responses to temperature also varied spatially (parameters β; for posterior mean estimates, see electronic supplementary material, table SC). The magnitude and direction of such responses were species specific—leaf out varied between 6.3 and 2.2 days earlier per degree Celsius (F. koreana), 4 and 2.6 days earlier per degree Celsius (G. biloba), 4.1 and 3.2 days earlier per degree Celsius (M. bombycis) and 7 days earlier and 9 days later (R. mucronulatum). First flowering also showed a range of responses, from 5.5 days earlier per degree Celsius to 6.4 days later (F. koreana), from 18.8 to 1.9 days earlier (P. persica), 3.8 days earlier to 5.8 days later per degree Celsius (P. yedoensis), 5.4 days earlier to 11.6 days later per degree Celsius (R. kaempferi), 4.5 days earlier to 10.4 days later per degree Celsius (R. mucronulatum), 10.6 days earlier to 18.4 days later per degree Celsius (V. mandshurica) and 2.2 to 1.6 days earlier per degree Celsius (Z. japonica). Such variability was also reflected in the predictions (figure 3). Species with more consistent patterns, e.g. M. bombycis, had tighter predictions than did species with large disparities between sites, e.g. P. persica. Exploration of the temporal random effects (how much of the year to year variability is explained by time; not shown) also differed among species. The slope of the regression time random effects versus year varied from downward trends (lower random effects with time) for F. koreana leaf-out (slope = −4.20 (s.e.: 0.45) and R2 = 0.5), with non-significant trends for most species, to upward trends (slope = 0.18 (s.e.: 0.02) and R2 = 0.6) for Z. japonica first flowering time. Plots of predicted versus observed phenology were also used to evaluate the fit of the model (electronic supplementary material, figure SA). These models showed a high-degree predicted fit of the observed values.

View this table:
Table 4.

Results: species level parameter estimates, values are posterior means ± s.d. (95% credible intervals) for model 6. Bold values indicate b coefficients that do not include zero in their credible intervals and denote a significant effect of temperature on phenology.

(ii) Autumn phenology

In the case of autumn phenologies, the species-level parameters associated with temperature, b, reveal a positive response to warmer temperatures—that is, they occur later in warmer years (table 4 and figure 4). These values are greater in absolute value than those seen for spring phenology (table 4), and here again, we found a variety of responses among the different sites (electronic supplementary material, table SC). Changes in the timing of leaf colour change varied from 7.9 to 2.4 days later per degree Celsius (A. palmatum) and from 11 to 1 days later per degree Celsius (G. biloba). Leaf fall also varied along the sampled range for each species, going from 4.3 to 2.9 days later per degree Celsius (G. biloba), 12.8 to 1.69 days later per degree Celsius (Morus) and 14.2 to 1.49 days later per degree Celsius (P. yedoensis). Such variability also affected the species' overall predictions as it did with spring phenology (figure 4). Trends in time random effects (regression of time random effect versus year) were positive; that is they increased along the time series. The slope parameters ranged from 0.018 (s.e.: 0.015; R2 = 0.02) for G. biloba leaf colour change to 0.18 (s.e.: 0.01; R2 = 0.76) for A. palmatum. The plots of predicted versus observed phenology show again a high degree of predicted fit of the observed data (electronic supplementary material, figure SB).

Figure 4.

Autumn phenology records (circles) and predicted overall species autumn phenology response to temperature (black lines, mean (parameters a and b) and 95% predictive interval). Colours represent data (dots) and predicted local responses (lines (parameters αs and βs)) at two locations.

(iii) Forecasting phenology

We used the parameters shown above to estimate species' overall phenological responses to temperature and to forecast potential changes in phenology under global warming for particular species at specific locations (see figure 5 for an example). We used the parameters calculated at the species level (a, Embedded Image, b, Embedded Image, timemean, Embedded Image and Embedded Image) to estimate the species' response to the range of temperatures for which we have data (figures 3 and 5, black lines). The forecasts represent range-wide estimates of the species' likely phenological response to the gradient of temperature values. The range of the response (95% credible interval) represents the natural variability in phenology found in the data, including unexplained spatial variation. Using the site-specific parameters (αs and βs) estimated through the hierarchical structure, we can now realistically extrapolate site-specific forecasts beyond the range of temperatures collected at that site (figure 5, red lines). These site-specific predictions greatly differ from those calculated using a similar linear model and the site's data in isolation (figure 5, blue lines).

Figure 5.

Comparison of forecasted phenological responses to a temperature gradient using a hierarchical approach and a non-hierarchical approach. Black lines indicate the overall species response (mean and 95% credible interval) across all sites from the hierarchical analysis. Black dots indicate observed phenology for a particular site. Red lines indicate the predictions (mean and 95% CI) for that site using the hierarchical approach. Blue lines indicate predictions (mean and 95% confidence interval) using a non-hierarchical approach, i.e. simple linear regression based only on data from that site, a technique commonly used in other studies. Vertical lines show the average spring temperature at the site (solid line) and predictions under two climate scenarios (dashed lines and shaded area).

4. Discussion

The combination of the hierarchical approach and the temporal and spatially extensive dataset allowed us to overcome the two major limitations that arise in the analysis of historical phenology data. First, using data on temporal and spatial variability in phenology, we were able to estimate each species' overall response to temperature while accounting for unexplained spatial variability. This allowed more reasonable forecasts for areas lacking data. Second, we used the large spatial extent of data in this study to more realistically forecast phenological trends at particular sites beyond their observed range of temperatures. For example, trends in southern locations did not determine our predictions in colder sites, but helped to inform the overall response under warmer climate scenarios.

Phenological responses to temperature may not have been evident at particular sites, based on simple linear regressions restricted to data at that site (figure 5, blue lines). This result, however, does not necessarily mean that the phenology of the species in question is independent of temperature. Site-specific factors, including precipitation, shading, soil conditions, nutrient concentrations and pathogens, can act together with temperature to determine a species' phenology (Winder & Cloern 2010). Hierarchical Bayes allowed us to account for this complexity through the use of site effects and latent (unobserved) variables. Modelling site effects hierarchically maintained variability that arises from unknown differences among sites, preventing asymptotic decline of variability owing to sample size (Clark 2005). The inclusion of latent variables allowed us to admit uncertainty caused by the difficulty of observing those variables. For example, including temperature as a latent variable helped account for differences between weather station measurements and temperatures actually experienced by individual plants.

A meta-analysis of European phenological events (Menzel et al. 2006) reported an overall advancement in spring phenology of 2.5 days for each degree Celsius, and a delay of autumn phenology of 1 day for each degree Celsius. In our dataset, for each degree Celsius increase in temperature, spring phenological events occurred 0–8 days earlier on average and autumn phenological events occurred 4–5 days later. The advances in springtime events and delays in autumn events were broadly consistent across species and locations in our dataset. However, despite this broad pattern, some species exhibited substantial spatial variability.

Most studies of long-term changes in phenology have highlighted advances in the timing of spring events and underemphasized the variability in those responses. However, some studies have shown trends towards later spring phenology in plant and animal species in Japan and other locations (Fitter et al. 1995; Butler 2003; Doi 2008). Such contrasting results are likely to be found more frequently as investigators continue to explore spatial variation in phenological trends. Here, it appears that spatial variability was probably driven not by latitude, but rather by other site-specific factors such as precipitation, soils, biotic interactions, human development, etc. (de Beurs & Henebry 2003; Touchon et al. 2006; Franks et al. 2007). Later spring phenologies in some locations may also be explained in part by species' inability to meet chilling requirements under warmer conditions, thereby delaying flower or leaf production (e.g. Morin et al. 2009).

The overall advances in spring plant phenology agree with findings in other plant species in East Asia (Kai et al. 1996; Chen 2003; Doi & Katano 2008), but contrast with the general lack of changes in spring animal phenology in the region (Doi 2008; Doi et al. 2008; Primack et al. 2009). The striking differences between phenological changes in plants and animals in the region are potentially reason for concern. It is possible that key plant–animal interactions, such as pollination, herbivory and also animal–habitat interactions (e.g. bare ground for nesting birds after snowpack melt), could be disrupted by rapid changes in phenology (Stenseth & Mysterud 2002; Visser & Both 2005). Evidence for such temporal mismatches has been found in other, more intensively studied locations (Inouye et al. 2000; Both et al. 2006; Post & Forchhammer 2008), but not in East Asia. Additionally, the spatial variation in phenological changes that we have observed could contribute to the potential for temporal mismatches, particularly for migratory species or species with large ranges (Post et al. 2008). We believe that the potential for phenological mismatches is an area deserving of much new research in East Asia.

As a result of earlier spring phenology and later autumn phenology, the growing season is lengthening in East Asia, as it is elsewhere in the world (Keeling et al. 1996; Menzel & Fabian 1999; Gordo & Sanz 2009). However, the relative importance of changes in spring and autumn phenology may vary around the world. In Japan and South Korea changes in autumn phenology appear to have a greater effect on growing season length than changes in spring phenology, in contrast to observations from Europe (Menzel & Fabian 1999; Menzel et al. 2006; Gordo & Sanz 2009). This stronger response of phenology in the autumn compared with spring was also seen in an earlier study of Ginkgo in Japan (Matsumoto et al. 2003). It is not clear what factors contribute to this difference between phenological changes in Europe and East Asia. The European studies (Menzel et al. 2006; Gordo & Sanz 2009) of changes in autumn phenology comprise a variety of ecosystems ranging from boreal forests to shrub semi-deserts. It is possible that the importance of precipitation for phenology in arid European sites contributes to the difference. Regardless of the cause, the disparities between advanced spring phenology and delayed autumn phenology could lead to very different responses in community dynamics and ecosystem processes. This topic deserves further investigation.

In many cases, investigators have used remotely sensed data from satellites to provide phenological data over broad areas (Myneni et al. 1997; Zhang et al. 2004). However, because phenological responses to global warming are idiosyncratic across species and regions, and because it is unknown exactly what phenological events are captured by various indexes derived from remotely sensed images (Schwartz 1998; White et al. 2009), remotely sensed data are generally of limited value for making species-specific forecasts. Also, although temperature is a major driver of many phenological events (Cleland et al. 2007), it interacts with many other local factors to determine the physiological responses of the plants and their phenology. This complexity makes generalizations from one place to another very difficult, and highlights the importance of geographically extensive datasets and modern analytical approaches for discerning species' responses to a changing environment (Primack et al. 2009). If one's goal is to identify ecological mismatches that might endanger the persistence of certain species, we suggest that data be collected at the level of the individual, in the context of the particular ecosystem, and along the distributional range of the species. This collection scheme would allow analytical models to consider the range of temperatures and site effects that occur throughout the range of the species.

Using a spatially extensive dataset, we have shown here that hierarchical Bayes in particular will be useful for incorporating the complexity of species responses into forecasts for data-poor sites and under warmer conditions. Though, as discussed earlier, these forecasts still have limitations—e.g. it is difficult to make forecasts for areas without existing data and for novel climate scenarios (beyond the range of variability described in the dataset), the information these forecasts provide will be critical to anticipate the effects of future phenological trends under global warming.

Overall, warming temperatures are related to earlier spring and later autumn phenology, which in turn are leading to longer growing seasons. Nonetheless, some sites show the opposite patterns, and locations with similar trends still vary. This indicates that the phenological response to temperature is complex and interacts with other local conditions. Only a clear understanding of the mechanisms determining phenology together with an inferential approach that can take into account the data's variability will allow us to formulate reliable predictions of future trends. Such understanding will require not only spatially extensive datasets but also comprehensive information of the environment in which these species live, and as researchers we should aim our efforts to gathering such data. Further, these results from meteorological observation sites should be extended to natural forests and other ecosystems in East Asia. We will need to sample the variability inherent in natural populations if we want to reliably quantify the scope of the phenological response to global warming.

Acknowledgements

We thank J. Diez, D. Katz and Z. Brym and three anonymous reviewers for comments on the manuscript. Support for this project was provided by the US National Science Foundation (DEB 0842465) to R.P., J.S. and I.I., and by the National Research Foundation of South Korea (2009-1419-1-6) to S.L.

Footnotes

    References

    View Abstract