A 250-year index of first flowering dates and its response to temperature changes

Tatsuya Amano, Richard J. Smithers, Tim H. Sparks, William J. Sutherland

Abstract

Widespread concerns about global biodiversity loss have led to a growing demand for indices of biodiversity status. Today, climate change is among the most serious threats to global biodiversity. Although many studies have revealed phenological responses to climate change, no long-term community-level indices have been developed. We derived a 250-year index of first flowering dates for 405 plant species in the UK for assessing the impact of climate change on plant communities. The estimated community-level index in the most recent 25 years was 2.2–12.7 days earlier than any other consecutive 25-year period since 1760. The index was closely correlated with February–April mean Central England Temperature, with flowering 5.0 days earlier for every 1°C increase in temperature. The index was relatively sensitive to the number of species, not records per species, included in the model. Our results demonstrate how multi-species, multiple-site phenological events can be integrated to obtain indices showing trends for each species and across species. This index should play an important role in monitoring the impact of climate change on biodiversity. Furthermore, this approach can be extended to incorporate data from other taxa and countries for evaluating cross-taxa and cross-country phenological responses to climate change.

1. Introduction

Human activities in the modern age have caused a serious loss of global biodiversity. Consequently, the Convention on Biological Diversity has adopted a target of ‘achieving by 2010 a significant reduction of the current rate of biodiversity loss at the global, regional and national level’. Scientists and policy-makers need to establish indicators of biodiversity and ecosystem services to assess progress towards the 2010 target and to effectively monitor the status of biodiversity (Balmford et al. 2005; Walpole et al. 2009).

In recent decades, climate change has been recognized increasingly as one of the most influential drivers of changes in biodiversity, particularly the distribution and abundance of plants and animals, and the timing of biological events (Walther et al. 2002; Parmesan & Yohe 2003; Root et al. 2003; Thomas et al. 2004; Parmesan 2006). With concerns about climate change continuing to escalate (Intergovernmental Panel on Climate Change 2007), there is an urgent need to develop indicators that show the potential impact on biodiversity and ecosystem services. This is considered one of the highest priorities in relation to the 2010 target (UNEP-WCMC 2009). However, to date, with reference to the 2010 target, only one study has established an indicator of the impact of climate change on bird populations (Gregory et al. 2009), and no existing indicator summarizes community-level changes in the timing of biological events (Walpole et al. 2009). This situation is critical because changes in phenological events across trophic levels would cause the temporal disruption of species interactions, such as phenological mismatches between plants and their pollinators (Hegland et al. 2009) or between birds and their plant and insect food supplies (Visser & Both 2005), possibly leading to an increase in extinction risks and loss of ecosystem services (Both et al. 2006; Memmott et al. 2007).

Although community-level indicators of phenological changes would play an important role in summarizing, and aiding policy-makers and the public to visualize, the impact of climate change on biodiversity, developing such indicators is challenging for two reasons. Firstly, previous studies have relied on long-term records collected continuously from a single location to investigate changes in phenology (e.g. Fitter et al. 1995; Fitter & Fitter 2002; Miller-Rushing & Primack 2008). Despite efforts to collect a comprehensive set of data (Parmesan & Yohe 2003; Root et al. 2003; Menzel et al. 2006), long-term records are still rare and unavailable for many regions and species (Primack et al. 2004; Miller-Rushing et al. 2006). Secondly, many studies have documented large differences in phenological responses to climate change among species (Menzel et al. 2006; Parmesan 2007), which are difficult to summarize at a community level (Parmesan 2007; Primack et al. 2009).

This study applies hierarchical models to estimate indices of phenological change from existing data, including short-term records from many locations. Statistical approaches have demonstrated that these hierarchical models can successfully estimate indices of population sizes from data with many missing values (Link & Sauer 2002). Recently, hierarchical models have been increasingly adopted to obtain a community-level overview of species-level phenomena, while accounting for variations in responses among species (Dorrough & Scroggie 2008; Royle & Dorazio 2008; Zipkin et al. 2009). Hierarchical models assume that species-level parameters are random effects governed by hyperparameters, which are used as community-level indicators of dynamics of interest (Zipkin et al. 2009). This approach can be applied to estimating indices of phenological change while accounting for the lack of long-term continuous records and variations in trends among species.

This study focuses on first flowering dates in plants, which are known to: (i) have a wide range of ecological, agricultural and socioeconomic consequences (Penuelas & Filella 2001), (ii) have been affected by recent climate changes (Fitter & Fitter 2002; Miller-Rushing & Primack 2008), and (iii) have been suggested as an indicator of the impact of climate change (Menzel et al. 2006). This study first uses 395 466 observation records of annual first flowering dates for 405 species from 1753 to date at sites throughout the UK and applies a hierarchical model to estimate a 250-year index showing a historical trend in the time of first flowering. Then, using the records between 1891 and 1947, where records are in relative abundance for many species, this study investigates the sensitivity of the index to the number of species and records per species included in the model.

2. Material and methods

(a) Data

Data were supplied by the UK Phenology Network (www.naturescalendar.org.uk). A national network was previously run by the Royal Meteorological Society (RMS) from 1875 to 1947 and re-established by the Centre for Ecology and Hydrology (CEH) in 1998. The Woodland Trust joined forces with CEH in 2000 to bring the project to a far wider audience through the web. Today, approximately 40 000 people across the UK are registered as recorders, providing extensive geographical coverage. More than three million phenological records have now been entered onto the database, including the RMS dataset and a substantial number of longer-term datasets dating back to the eighteenth century (data source in appendix S1, electronic supplementary material). The data used in this study consist of 395 466 records of first flowering dates for 405 species from 1753 to date. Dates were converted to days of the year (1 = 1 January etc.) prior to analysis. If more than one record of a particular species occurred at the same location in the same year, the mean of all the records was used. To reduce calculation time, the maximum number of records per species per year was set at 300 and records were randomly selected, as necessary. Species included in the analysis, their observation periods and the total numbers of records and sites are listed in table S1, electronic supplementary material. The temporal changes in the numbers of species and records per species involved in the analysis are shown in figure S1, electronic supplementary material.

The Central England Temperature (CET) series (Manley 1974; Parker et al. 1992) has been shown to be broadly representative of temperatures in other parts of the UK (Croxton et al. 2006). The CET was downloaded from the webpage of the Meteorological Office Hadley Centre (http://hadobs.metoffice.com/hadcet/). The monthly mean CET was used to calculate the mean for each period.

(b) Model

The number of sites contributing to the flowering data, the observation period and temporal trends in first flowering dates can vary greatly among species (table S1, electronic supplementary material). Thus, the model estimates the community-level trend by accounting for: (i) the absolute difference in first flowering dates between species, (ii) the variation in the temporal trend of first flowering dates between species, and (iii) the latitudinal difference in first flowering dates within species.

First, the species-level index βit for species i in year t is assumed to be derived from a normal distribution with the mean of the community-level index βt showing a temporal trend in the first flowering dates shared by all the species as follows: Embedded Image 2.1 Embedded Image is a hyperparameter that allows the model to account for variations in the temporal trend of first flowering dates among species, the importance of which has been indicated by earlier studies (Parmesan 2007; Miller-Rushing & Primack 2008). The advantage of this hierarchical structure is that (i) the hyperparameter (βt in this case) can be used to summarize community-level phenomena (Dorrough & Scroggie 2008; Zipkin et al. 2009) and (ii) parameter estimates from individual datasets (βit in this case) can be improved over those that would be obtained by fitting separate models to each dataset (Jonsen et al. 2003; Zipkin et al. 2009). Second, the expected first flowering date μitj for species i in year t at site j is expressed using the species-level index, the species effect αi and the northing of the survey sites Nj as a covariate: Embedded Image 2.2 αi accounts for the absolute difference in first flowering dates among species, and is drawn from a mean zero normal distribution with variance Embedded Image, which is a hyperparameter. Here the use of normal distributions is supported by the symmetric shape of the histogram for mean first flowering dates of different species (figure S2, electronic supplementary material). bi is the species-specific coefficient of the effect of the northing, and is drawn from a normal distribution with mean m. b and variance Embedded Image, which are also hyperparameters. This structure enables the latitudinal difference in first flowering dates within species to be modelled. For simplicity, the difference in the response to temperature across latitudes was not considered because with some species, the number of sites was insufficient for estimating spatial variations in phenological responses. But spatial variability in phenological responses has been reported by earlier studies (Parmesan 2007; Primack et al. 2009), and thus, incorporating such effects should be an important challenge for the future. Finally, the observed first flowering date yitj for species i in year t at site j is assumed to follow a normal distribution: Embedded Image 2.3 Hierarchical models are analysed using Bayesian methods, which naturally account for uncertainty in parameter estimation (Royle & Dorazio 2008), and the model was fitted to the data with WinBUGS 1.4.3 (Spiegelhalter et al. 2003; script provided in appendix S2, electronic supplementary material) and a package R2WinBUGS in R (Sturtz et al. 2005; script in appendix S3, electronic supplementary material). Gamma distributions with mean of 1 and variance of 100 were used as prior distributions for the inverses of Embedded Image and Embedded Image. Normal distributions with mean of 0 and variances of 100 000 and 1000 were used as prior distributions for βt and m. b. Each MCMC algorithm was run with three chains with different initial values for 190 000 iterations, with the first 30 000 discarded as burn-in and the remainder thinned to one in every 40 iterations to save storage space. Model convergence was checked with R-hat values (Gelman et al. 2003) and trace plots of all the chains for sampling (Spiegelhalter et al. 2003).

(c) Analyses

The relationship between the index and CET was assessed using a simple linear regression. In a preliminary analysis, there was a negative relationship between residuals from the ordinary least-squares estimates and the number of (i) species (spt) and (ii) sites per species (sitet) used to estimate the index each year, which indicated non-homogeneity of the variance (Quinn & Keough 2002). Thus, a linear regression with weighted least squares was applied (Quinn & Keough 2002). To consider the effects of the number of (i) species and (ii) sites, the weight wt was defined as follows: Embedded Image 2.4

A linear regression model was fitted by minimizing Embedded Image, where Embedded Image represents the index value predicted by the regression line (Quinn & Keough 2002).

(d) Simulation tests

Since the number of species and records per species included in the model vary greatly throughout the period considered, simulations were conducted to test the sensitivity of the estimated index to the number of species and records per species included in the model. The hierarchical model was first applied to records between 1891 and 1947, when relatively many records are available on many species (figure S1, electronic supplementary material), to derive the community-level index, which was assumed as the ‘true’ index to be estimated in the following two simulations. First, 75, 50, 25 and 5 per cent of the species observed between 1891 and 1947 were chosen randomly, and the community-level indices were estimated using these subsets of species. Second, 75, 50, 25 and 5 per cent of the records for each species between 1891 and 1947 were chosen randomly and used for estimating the community-level indices. To explore the effects of the spatial distribution of records, indices were also estimated using data with 75, 50, 25 and 5 per cent of available records chosen by ascending order of northing. Each procedure was repeated 10 times and correlation coefficients between the true and estimated indices were calculated.

3. Results

Figure 1 shows the 250-year community-level index of the time of first flowering from the applied hierarchical model. There is a clear advance in the time of first flowering in recent decades. The mean community-level index for the most recent 25 years indicated an earlier flowering than the mean for any other consecutive 25-year period between 1760 and 1984, and the difference ranged from 2.2 (95% credible interval: 1.0–3.4; versus 1910–1934) to 12.7 days (9.8–15.6; versus 1835–1859, figure 1). Although the index showed a close negative correlation with the December–February, January–March and March–May means of CET (figure S3, electronic supplementary material), the correlation was strongest for the February–April mean (figure 2, weighted regression coefficient =− 4.96, p < 0.001, R2 = 0.58). The rate of advance in the index with temperature (5.0 days earlier per 1°C) corresponds approximately to the range reported for a wide variety of species (Sparks et al. 2000; Miller-Rushing & Primack 2008). The record for the most recent 25 years clearly lies in the area of high temperature and early-first flowering (figure 2). Figure 3 shows the index and February–April temperature for the last 250 years. The model also estimated species-level indices, which showed, for example, that the time of first flowering for both hawthorn Crataegus monogyna and blackthorn Prunus spinosa occurred earlier in the most recent 25 years than any other consecutive 25-year period during the 250- and 125-year observation periods, respectively (figure 4).

Figure 1.

The median (red line) and 95% credible intervals (grey area) of the estimated community-level index (day of the year) showing a temporal change in the timing of first flowering shared by 405 plant species observed throughout the UK. The black line indicates the mean for every 25 years and the dotted line that for the most recent 25 years. The years without estimates indicate those without any observation records (1766, 1813, 1814 and 1817).

Figure 2.

The significant negative correlation between the community-level index and the mean of the CET for Feb–Apr. The records for the most recent 25 years (red points) lie in the area of high temperature and early-first flowering.

Figure 3.

The community-level index (red) and the mean of the CET (blue) for Feb–Apr (index inverted for clarity).

Figure 4.

Example of the species-level indices showing changes in first flowering dates of (a) hawthorn and (b) blackthorn. Lines and areas are as defined in figure 1. In this figure, to show the absolute difference in first flowering dates among species, the species effects (αi in equation (2.2)) are added to the species-level indices. The years without estimates indicate those without any observation records (hawthorn: 1766, 1811–1817, 1873, 1883, 1948–1951 and 1960–1973; blackthorn: 1753–1808, 1812–1863, 1870–1890 and 1945–1951).

The simulation tests using records from between 1891 and 1947 revealed that the indices estimated with the reduced number of species were closely correlated with the true index, which was estimated with all the records, even when just 25 per cent of all the species were included (r > 0.8; figure S4, electronic supplementary material). However, with just 5 per cent of the species, the correlation coefficients between the estimated and true indices varied greatly, from almost zero to over 0.8 among different subsets of species (figure S4, electronic supplementary material). On the other hand, the proportion of records per species included in the model did not seem to have a great impact on the estimates; the indices estimated even with 5 per cent of the original records per species were closely correlated with the true index (r > 0.8), and the result was the same when the records were chosen in ascending order of northing (figure S5, electronic supplementary material).

4. Discussion

The results of this study have implications for both science and policy.

Firstly, the hierarchical model adopted in this study allows species-level observation records from more than one site to be integrated into a single community-level index, as well as species-level indices. Earlier studies have relied on data collected continuously from a single location to investigate phenological responses to climate change (Fitter et al. 1995; Fitter & Fitter 2002; Miller-Rushing & Primack 2008). The limited number of long-term records of phenological events has restricted research on the topic in most areas of the world and for many species (Miller-Rushing et al. 2006). However, the model developed in this study makes the best use of existing data, including short-term records from many locations, to create an index. This is similar to the approach used by Kauserud et al. (2008) to derive a community-level phenological response of mushrooms to climate. A critical difference is that it accounts for variations in temporal trends among species, which if ignored can lead to inappropriate conclusions about community-level trends. This structure also enables species-level indices to be estimated. Thus, the approach developed in this study should promote understanding of the impact of climate change on ecological systems for previously unexamined localities and species.

Further, given the gap in temporal/spatial scales between ecological and climatological studies (Harrington et al. 1999; Walther et al. 2002), another advantage of this study is that it has estimated an index of first flowering date over a 250-year period covering the whole country. Consequently, for example, this study revealed that the mean index of first flowering for 1910–1934 was close to, though not earlier than, the mean for the most recent 25 years. Since conclusions about the impact of climate change based on short-term observations can be misleading (Chuine et al. 2004; Kerr 2009), the 250-year index developed in this study should provide an important context for investigating ongoing phenological responses to climate change.

Although satellite observations and climate-based models are becoming increasingly important for studies of phenological changes over space and time (Schwartz et al. 2006; Cleland et al. 2007; Morisette et al. 2009), one major barrier to progress in such indirect approaches is a lack of appropriate large-scale ground phenological data that could be effectively compared with metrics derived from satellite/climate data (Cleland et al. 2007). Our study has successfully scaled up ground phenological data from site-level to national-level indices and thus, should serve as a benchmark for the comparison with outputs from other approaches.

The simulation tests revealed that the estimated community-level index is relatively sensitive to the number of species, not records per species, included in the model. Figure S4, electronic supplementary material, indicates that the index value based on only a small proportion of the species in the target community is strongly dependent on the composition of species used for the estimation. Records from more than 5 per cent of all the species, chosen randomly, are required for the index to reflect the temporal trend of the whole community. Thus, for instance, special care is needed when interpreting the community-level index for the eighteenth and nineteenth centuries because for most of this period, the index reflects the temporal trend of only a few species included in the model (figure S1, electronic supplementary material). Nevertheless, from the late nineteenth century onwards, relatively many species were involved in the estimation of the community-level index. New records can also be added to the analysis as new datasets become available. This may become particularly important when additional records are discovered for the eighteenth and nineteenth centuries. A key investigation for the future will be to find out which species show similar (and dissimilar) trends with the community-level index. Identifying such species would help (i) avoid unintentional reliance only on some species with exceptional trends and (ii) determine the combination of species that allows the effective community-level index to be estimated with minimum effort. Moreover, while this study was based on first flowering dates as the only available data source, which can be affected by both population size and sampling frequency (Miller-Rushing et al. 2008), future work should carefully evaluate the effects of these factors on consequent indices.

Secondly, this study has successfully developed an index that shows a community-level phenological response to climate change that can be adopted by governments for monitoring the impact of climate change. Although there is now a broad consensus that climate change poses a threat to global biodiversity, efforts to develop indices that indicate the impact of climate change on ecosystems have been scarce (Gregory et al. 2009; Walpole et al. 2009), creating an urgent need to develop an index with reference to the 2010 target (UNEP-WCMC 2009). The importance of developing community-level indicators lies in summarizing, and aiding policy-makers and the public to visualize, the impact of climate change on biodiversity, as is also achieved by the indicators of global average temperature, sea level and Northern Hemisphere snow cover (Intergovernmental Panel on Climate Change 2007).

At the moment, there is only one phenological index that has been officially adopted by a government: the UK Spring Index adopted by the British government (Department for Environment Food and Rural Affairs 2009). However, this index simply calculates the annual mean observation dates of just four biological events: first flowering of hawthorn, first flowering of horse chestnut Aesculus hippocastanum, first recorded flight of orange-tip butterfly Anthocharis cardamines and first arrival of swallow Hirundo rustica (Department for Environment Food and Rural Affairs 2009). The approach adopted by the UK Spring Index requires a complete set of long-term records for each species, which limits the range of applicable species. This seems to be a critical problem for summarizing the long-term community response because (i) the index, based on only a few species, may not accurately reflect the trend of the community, as shown in our simulation and (ii) even a single missing value for one species in a given year would make it impossible to calculate the index for that year. In contrast, this study is based on phenological trends for 405 species observed throughout the country and thus, has successfully developed a long-term index that summarizes the community response.

Given the presence of numerous phenological observation networks all over the world (Cleland et al. 2007), future efforts should be directed towards creating indices across taxa and regions. Our new model can be easily extended to incorporate records of other taxa or those in other regions of the world. Further, to monitor the level of disruption in species interactions among trophic levels—one of the processes through which climate change affects biodiversity (Visser & Both 2005)—it will also be important to develop and compare indices for different trophic levels. Therefore, the approach developed in this study provides a novel and powerful tool for monitoring the impact of climate change on ecological systems at global, regional and national levels.

Acknowledgements

We thank Y. Yamaura for comments on early drafts, S. V. Prior for preparing tables, more than 100 000 volunteers for collecting the data and M. Amano for all the help. T.A. is funded by Grant-in-Aid for Young Scientists B (21710246), Japan Society for the Promotion of Science and W.J.S. by Arcadia. We thank referees for their comments on an earlier version of this paper.

Footnotes

    • Received February 11, 2010.
    • Accepted March 16, 2010.

References

View Abstract