Climate change affects the phenology of many species. As temperature and precipitation are thought to control autumn color change in temperate deciduous trees, it is possible that climate change might also affect the phenology of autumn colors. Using long-term data for eight tree species in a New England hardwood forest, we show that the timing and cumulative amount of autumn color are correlated with variation in temperature and precipitation at specific times of the year. A phenological model driven by accumulated cold degree-days and photoperiod reproduces most of the interspecific and interannual variability in the timing of autumn colors. We use this process-oriented model to predict changes in the phenology of autumn colors to 2099, showing that, while responses vary among species, climate change under standard IPCC projections will lead to an overall increase in the amount of autumn colors for most species.

Funding: Research at Harvard Forest is partially supported by the National Foundation′s (NSF) LTER program (award number DEB-0080592). ADR acknowledges support from the National Science Foundation through the Macrosystems Biology program (award EF-1065029), and the Northeastern States Research Cooperative, in addition to the United States Geological Survey Status (USGS) and Trends Program, the United States National Park Service Inventory and Monitoring Program, and the United States of America National Phenology Network, through grant number G10AP00129 from the USGS. The contents of this paper are solely the responsibility of the authors and do not necessarily represent the official views of either USGS or NSF. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Copyright: © 2013 Archetti et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

In order to predict how autumn colors may respond to forecast changes in environmental drivers, we analyzed data on leaf color change collected annually between 1993 and 2010 in a New England forest for eight study-species that develop anthocyanins in autumn. For each species we calculated the average percentage of colored leaves and of fallen leaves for each day of the year for the 18 years during which the data were gathered. We investigated correlations between temperature and precipitation during different times of the year, and the timing of various autumn color thresholds and leaf fall dates. We compared two types of models to explain autumn coloration and leaf fall. First, we used an empirical approach [26] based on stepwise multiple linear regression, with monthly means of temperature and precipitation as the candidate independent variables. Second, we used a more mechanistic approach using a cold-degree-day photoperiod-dependent model [27] . The correlation analysis and empirical modeling allow us to identify environmental drivers that may be missing from the mechanistic model, which is highly constrained in its structure, and which does not, for example, account for relationships between precipitation and autumn color. We evaluated the models against the observational data using cross-validation methods. We then used the most robust modeling approach, in conjunction with IPCC climate projections, to forecast changes in the phenology of autumn color and leaf fall, between now and the year 2099.

Additionally, the potential impact of climate change on the intensity and duration of autumn coloration is, in some regions, of enormous economic importance [25] . Autumn tourism—much of which is to participate in so-called ‘leaf peeping’—contributes billions of dollars each year to the economies of the states of the eastern U.S.A. and provinces in adjacent Canada. If climate change reduces the duration of autumn color display, or results in less vibrant displays, future tourism revenues will likely be reduced.

At the continental scale, warmer autumns have for instance been related to lower net carbon fixation [20] – [21] , as a consequence of a higher enhancement of ecosystem respiration than the concomitant enhancement of gross photosynthesis. At a local scale, temperate deciduous forests may on the contrary show a higher annual net carbon fixation during warmer autumn as a consequence of an extended leafy season [22] . There is further evidence that the asynchrony of autumn phenology may alter the competition between co-occurring plant species, either in the case of symmetric (between understory plants - all plants being light-limited by the overstorey canopy) [23] or asymmetric (between overstory and understory plants) [24] competition.

About 15% of the tree species of the temperate regions of the world change their leaf color from green to yellow or red in autumn, a percentage that can reach 70% in some regions like New England (Northeast USA) [14] – [15] . As leaf color change and leaf fall are thought to be controlled by temperature and precipitation [16] – [18] , it is possible that climate change may also affect autumn phenology, with obvious biological and ecological implications [19] .

Temperature affects biological processes ranging from the molecular to the ecological level. It is not surprising, therefore, that climate change is altering the phenology of many species [1] – [7] . In plants, the impacts of climate change on spring phenology (flowering) are well documented [8] – [13] . Much less is known, however, about how warming temperatures and altered precipitation regimes affect autumn phenology, specifically as related to leaf coloration and senescence.

Materials and Methods

Data We analyzed data on the autumn phenology of Acer rubrum (red maple), Acer saccharum (sugar maple), Fraxinus americana (white ash), Nyssa sylvatica (black gum), Prunus serotina (black cherry), Quercus alba (white oak), Quercus rubra (red oak) and Quercus velutina (black oak) at Harvard Forest, a research area owned and managed by Harvard University, in Petersham, Massachusetts, USA (Prospect Hill Tract; 42.54 °N, 72.18 °W). For more than twenty years, phenological observations have been made, every 3–7 days in spring and autumn [18], [28], by the same observer. The observed trees (3 to 5 permanently-tagged individuals per species) are located within 1.5 km of the Harvard Forest headquarters at elevations between 335 and 365 m above sea level. The field protocol for autumn observations was finalized in 1993 and here we use observations through the end of 2010. Beginning in September, and continuing through the end of leaf fall, leaf coloration (the percentage of leaves that have changed color on a given tree) and leaf fall (the percentage of leaves that have fallen from a given tree) are estimated for each individual observed. The raw data are available at http://harvardforest.fas.harvard.edu/data/archive.html (datasets HF000, HF001, HF003); the transformed data and the codes used for the analysis are available from the authors, while the final data are in Table S1.

Measures of autumn color We used the original data to infer the day (c x ) on which the percentage of colored leaves is x and the day (f x ) in which the percentage of fallen leaves is x (where x may take a value of 10, 25, 50, 75 or 90 percent). Assuming that both color and leaf retention change as a linear function between the days in which the observations were recorded, we derived c x using the formula where x INF and x SUP are the available measure immediately lower and higher than x; f x was derived in a similar way as For a few species, in some years (18 in a total of 2304 data points, that is 0.65% of the data), certain thresholds (mainly c 10 and c 25 ) had already been reached before the first field observations were made: in these cases, rather than extrapolate backwards, we simply treated these as missing data. We also used c x and f x to build two different measures of abundance of autumn color: d x = f 90 −c x measures the duration of autumn color as the number of days between the day when a percentage x of the leaves are red (c x ) and the day when 90% of the leaves have fallen (f 90 ). The amount of autumn color is measured by (i n −i n−1 )y n−1 +(i n −i n−1 )(y n −y n−1 )/2 if y n >y n−1 and by (i n −i n−1 )y n +(i n −i n−1 )(y n−1 −y n )/2 if y n <y n−1 , where y n = r n (1−t n /100); r n is the percentage of red leaves, t n is the percentage of leaves retained, i n is the (julian) day when the nth measure (of a total of m measures) was taken. The yearly amount of autumn color therefore is (in a Cartesian plane), the area below the lines that connect the daily amount of autumn color (see Figure 1). 100 units of A correspond to one calendar day in which all leaves are retained and red. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Figure 1. The amount of autumn colors over time for eight deciduous broadleaf species that turn red in autumn. The amount of autumn color (0–100) is calculated as i n (100−j n ) on day n, where the percentage of red leaves i n is multiplied by the percentage of leaves retained (100−j n ). Individual years (1993–2010) are shown by dotted lines, and their average by the thick curve. https://doi.org/10.1371/journal.pone.0057373.g001

Correlation analysis and regression modeling Air temperature and precipitation are measured (daily) at the Harvard Forest near to the trees on which phenological observations have been conducted. Data for the Shaler (1964–2002) and Fisher (2001-present) meteorological stations are available online at the web address given above; any missing observations were filled using measurements from the Harvard Forest EMS AmeriFlux tower, approximately 1 km distant. For both temperature and precipitation, we first calculated averages (of the daily measures) over all the 1- to 52-week timeframes preceding each day of the year. We then calculated the correlation coefficients between these averages and each of the measures of autumn color (A; c x , f x , d x ; see above) for each species. Based on the correlation analysis, we identified the periods of the year during which the largest positive and the largest negative correlations were observed with the measures of autumn color. For our empirical modeling of leaf color threshold dates (c x ) and leaf fall threshold dates (f x ), we calculated monthly means of temperature and precipitation during the leaf-on (May to October) period. We conducted a stepwise multiple linear regression procedure with the monthly mean drivers as candidate independent variables (6 months×2 drivers = 12 candidate variables). We specifically chose a monthly time interval (rather than weekly) for averaging, and restricted our analysis to the leaf-on period, so as to avoid having too many candidate variables, which could increase the likelihood of type 1 (false positive) errors and potentially lead to the inclusion of spuriously correlated variables in the regression. At each iteration of the stepwise procedure, variables that would be significant at a p-value of ≤ 0.20 were added to the regression but were subsequently removed if, after other variables were accounted for, the p-value exceeded 0.05. We fit a separate model to each c x threshold and each f x threshold; A and d x were then calculated from c x and f x . Below, we refer to this Multiple Linear Regressions approach as the MLR model.

Process-oriented modeling We used a cold-degree-day photoperiod-dependent (CDD/P) model [27]. This model was initially designed to simulate a coloring stage and was further applied in this study to the simulation of a fall stage. Whatever the senescence stage (c x or f x ) considered, it is defined in the model by S sen (arbitrary units) for each day (doy) following D start (the date at which a critical photoperiod P start is reached), representing the progress of the simulated process. Leaf coloring or fall reaches a given stage (c x or f x ) when S sen reaches a threshold value (Y crit , arbitrary units). In this model, the time derivative of the state of senescence (R sen , arbitrary units) on a daily basis is formulated as: Where P(doy) is the photoperiod expressed in hours on the day of year doy; T(doy), the daily mean temperature (°C); T b , the maximum temperature at which the considered senescence (i.e. coloration or fall) process is effective (°C); f[P(doy)], a photoperiod function that can be expressed as follows : The complete model therefore includes five parameters (P start , T b , x, y, Y crit ). The dummy parameters x and y may take any of the {0, 1, 2} discrete values, to allow for any absent/proportional/more than proportional effects of temperature and photoperiod to be included. A feature of this model structure is that, depending on the value of x, the modeled phenophase can be considered as dependent (x>0) or independent (x = 0) on cold-degree days. In the latter case, the occurrence of the phenophase is only determined by a threshold photoperiod. The optimization procedure consisted of exploring the whole space of parameters for P start (from 10 to 16 h with a 0.5 h step), T b (from +7 to +30°C with a 0.5°C step), x, and y. The Y crit parameter was identified through the Powell (gradient descent) optimization method [29]. Parameter optimization was based on minimizing the model-data mismatch, quantified in terms of root mean squared error. As with the MLR approach, the CDD/P model was fit independently on leaf color (c x ) and leaf fall (f x ) data for each species. Yet, while the MLR approach was fit on each color and fall stage (e.g. 5 fits for color from c 10 to c 90 ), we fit the CDD/P model over the complete phenological trajectory (e.g. simultaneously for all five stages from c 10 to c 90 for leaf coloration) defining for each model structure a set of five Y crit parameters, one per observed stage. We thereafter used the two CDD/P models fit independently on coloring and fall data to predict canopy duration (d x ) and the amount of color (A). Statistics were computed using MATLAB version 7.10 (The MathWorks Inc., 2010).

Robustness assessment of the modeling approaches The accumulation of a large phenological dataset requires sustained effort over many years, which is why multi-decadal records are relatively scarce. With 18 years worth of data, the Harvard Forest dataset is one of the longest autumn datasets published [28]. However, it is certainly possible that either the statistical (MLR) or process-oriented (CDD/P) approaches could result in models being over-fit to what is still a relatively short time series. After performing a first fit of both approaches on the full dataset, we evaluated the robustness of each model (i.e. the ability of the model to predict an unknown dataset) by using cross-validation analysis [30], [31]. This approach is commonly used when wholly independent data (e.g., from another site) are unavailable for model testing (for examples in the phenology literature, see [26]). Specifically, we used a one-out cross-validation, which is particularly appropriate when the dataset is relatively small. To conduct the cross-validation, the models were fit sequentially on 17 of 18 points (i.e. years) from the original dataset (‘calibration’) and tested for their ability to simulate the remaining point (‘validation’). This was repeated 18 times, so that each data point was included in the validation set exactly once. Model performance statistics (root mean square error, RMSE, and model efficiency, ME [32]) were then calculated across the 18 validation points. We assessed the ability of each of the two modeling approaches to maximize the trade-off between model parsimony and goodness-of-fit using Akaike's information criterion, corrected for small samples (AICc [33]).