Robust Model Predicts Shoot Phenology of Fraser Fir under Extreme Conditions

Fraser fir (Abies fraseri [Pursh] Poir.) is an important Christmas tree species in the United States, and understanding its phenology is important for managing Fraser fir trees in plantations or forests. Many management decisions are informed by and dependent on shoot phenology, from late spring frost protection to shearing, and from timing pesticide sprays to managing cone production. The ability to predict important phenological stages will become increasingly important as the climate warms, as is predicted for the primary regions where Fraser fir is grown for Christmas trees. Here, we report on the development of a model of shoot phenology in Fraser fir, and present one example of how this model may be applied to the problem of managing cone production. We surveyed shoot phenology at nine Christmas tree plantations in Michigan over three years, and used the data obtained to develop a phenology model of shoot growth. Derived from the beta sigmoid function and based on growing degree days, this phenology model offers a high predictive power and is robust to extremes of temperature and precipitation. When applied to cone production, our model provides guidance for timing practices that influence cone bud formation, both for reducing nuisance cones in Christmas tree plantations and for enhancing cone production in seed orchards. In addition, the model may assist with timing other practices tied to shoot phenology. The performance of our model under extreme heat and drought conditions suggests a role for this and other phenology models in predicting and mitigating the effects of climate change on tree growth and development.


Introduction
Fraser fir (Abies fraseri [Pursh] Poir.) is an important Christmas tree species in the United States, with annual sales of $100 million in the southern Appalachians [1], and a growing market share in the Midwest.Natural stands of Fraser fir are only found at high elevation sites in the southern Appalachian Mountains, where temperatures are cool and cloud-immersion is common [2].Fraser fir trees grown in Christmas tree plantations experience warmer and less humid conditions in the Appalachians, and more so in the Midwest [3].Climate models predict increasing temperatures and decreasing summer precipitation for both regions [4], raising questions about how growth and phenology will change in response, in addition to the effects this will have on the production and sale of this commercially important species [3,5].
Understanding phenology is important for managing trees grown in plantations or forests.For example, both balsam gall midge (Paradiplosis tumifex Gagne) [6] and balsam twig aphid (Mindarus abietinus Koch) [6][7][8] cause economically significant damage to Fraser fir trees, but insecticides must be properly timed, based on insect and tree phenology, to provide acceptable control.Similarly, many Christmas tree producers in the southern Appalachians apply the herbicide glyphosate for weed suppression between rows of trees, in a process termed 'chemical mowing' [9].The application rate and method vary during the growing season based on shoot phenology, to prevent injury to the expanding shoots and needles [10].In the United States, most growers shear their trees annually to produce the dense, conical form desired by consumers.The timing of shearing varies among producers, but is governed by phenological changes, such as terminal budset or cessation of lateral growth [11].Climate change is expected to usher in additional challenges [5].Insect pest density is predicted to rise with the warming temperatures [8].Early warming in the spring encourages a flush of tender growth that is sensitive to frost events later in the spring [12,13].
In Abies, strobilus bud initiation and differentiation are closely tied to shoot phenology, and occur in early summer, the year before cones emerge.Cone buds initiate from lateral buds when the growth of lateral shoots begins to slow.Differentiation of reproductive organs follows over the next two weeks, and is complete by the time lateral shoot expansion has ended [14][15][16][17][18]. Environmental signals during the period of initiation and differentiation interact with endogenous controls to regulate reproductive development [14].Gibberellin (GA) signaling is involved in the initiation of reproductive development [19], and GA is commonly applied to induce or enhance cone production in seed orchards [14,20,21].High temperatures and drought increase cone production in many forest trees [14,22], and techniques to increase tree stress, such as girdling or root pruning, are used in seed orchards to enhance cone production-generally in combination with GA [20,23].Tree response to GA, heat, and other treatments is dependent on tree phenology and therefore varies based on timing [24].
In contrast, heavy cone production in Fraser fir Christmas tree plantations is a significant problem for producers in the Midwest [25], and increasingly in North Carolina [26].Cones decrease the value of a tree by displacing lateral branches in the upper third of the tree crown, resulting in sparse tops that are less acceptable to consumers.Expanding cones also compete with vegetative growth for photosynthates [27,28].Therefore, many growers remove cones by hand while they are still small.Per tree, cone removal represents the highest labor expense for many growers [26].
Various experimental approaches have been undertaken to mitigate the effects of heavy cone production in Fraser fir Christmas tree plantations, all of which depend on shoot phenology.Cones may be mechanically dislodged using special tools, but this must be done prior to vegetative budbreak or new shoots will be damaged.Similarly, caustic sprays (chemical thinners or herbicides) may be used to abort cones in the spring while they are still small [25,26], but many will readily damage emerging foliage when applied at rates high enough to abort cones.An alternative approach is to disrupt strobilus initiation or differentiation.This approach may allow more lateral buds to develop into shoots, resulting in a more uniform branch density throughout the tree crown.Strobilus development may be disrupted by chemical treatments, such as the application of plant growth regulators (PGRs) [25,29], and cultural practices that modify environmental conditions [30].However, treatments need to precede reproductive bud initiation and continue throughout differentiation.PGR sprays, in particular, need to be carefully timed in order to be effective.For cone reduction, most relevant PGRs are GA inhibitors that interfere with gibberellin biosynthesis.Therefore, they might over-regulate stem elongation if applied too early, but may not affect cone production if applied too late.
Phenology is key to the adaptation of trees to their environment, and phenology models are important for predicting the impacts of climate change on the tree health, growth, and productivity [31].Shoot growth phenology advances when a plant is exposed to temperatures above a threshold base temperature that is required for that stage of development.Growing degree days (GDD) are the units used to represent this accumulation of temperature, or heat units, that govern growth and development over time.Therefore, shoot growth phenology can be better modeled using GDD rather than calendar date [32,33].In this paper, we report on the development of a GDD model of shoot phenology that may be helpful for timing insecticides, herbicides, PGR sprays, or cultural treatments for the management of Fraser fir.We then provide an example of how this model may be applied to specific problems related to cone production, by mapping the timing of reproductive bud initiation and differentiation, derived from the literature, to our model.The performance of our model under extreme heat and drought conditions suggests a role for this and other phenology models in predicting the effects of climate change on tree growth and development.Although developed in the Midwest, this model may be of interest to plantation and forest managers in the southern Appalachian Mountains, where the ability of Fraser fir to adapt to climate change is of economic and ecological import [3].

Materials and Methods
We divided the study into two phases.In Phase I, we developed a shoot phenology model based on terminal shoot growth measured at nine sites in Michigan (Figure 1).We selected leader growth for developing the phenology model because it can be easily and accurately measured with little specific training for grower-cooperators.In Phase II, we adapted the terminal shoot model to lateral shoot growth because lateral growth is closely tied to reproductive bud development and differentiation but lags behind terminal shoot growth.1.

Phase I: Development of Phenology Model
From May 2011-July 2012, we monitored terminal leader growth on Fraser fir trees in nine operational plantations representing a wide range of locations and site conditions in Michigan (Table 1).We randomly selected 25 trees at each site, and measured leader length each week beginning at budbreak and ending when the average leader length for all trees measured in a field was unchanged from the previous week.We paired this growth data with GDD data (base 5 • C, Baskerville-Emin method) obtained for each measurement date from nearby automated weather network stations operated by Michigan State University Enviro-weather (http://enviroweather.msu.edu).We developed our phenology model from the measurements of leader growth.Because the graph of the response variable (leader length) against GDD followed a sigmoidal curve, we initially fit a logistic function, which is commonly used to model biological growth data [36].Despite the high R 2 (0.95), we noted bias in the graph of residuals against predicted values, indicating that the model did not accurately fit the data at the beginning and end of the growth curve.We then fit a beta sigmoid function (BSF), which is a more flexible, generalized polynomial equation that is able to accurately represent a variety of sigmoid plant growth patterns [36][37][38], following the equation given by Shi et al. [37]: where t b ≤ t ≤ t e , l is shoot length as a percent of the maximum length, t is thermal time ( • Cd), and c m is the maximum growth rate occurring at time t m .This model is simplified to assume that the growth rate is 0 at the beginning (t b ) and end (t e ) of the growth period.We fit the model in R 3.2.4[39] using functions developed by Shi et al. [37].We trained the model on growth data from 2011 collected at five sites, located an average of 6.9 mi.(11.1 km) from Enviro-weather automated network stations that supplied the GDD data (Figure 1).Two data sets were used for external validation to test the predictive power of the model.The 2011 validation data set included growth data from 2011 collected at four study sites.For three sites, GDD data was calculated as an average of data from two or three Enviro-weather stations.For the fourth site, GDD data was obtained from a single weather station located nearby.The 2012 validation data set was primarily for temporal validation, and consisted of data from 2012, obtained from the five sites that contributed to the training data set and one site that contributed to the 2011 validation data set.
The phenology model was fitted and validated using the means of leader length data for each site and measurement date.To adjust for differences in vigor between sites, terminal leader length data were normalized as percent of maximum growth for each site by dividing the average leader length for each date by the average final leader length and multiplying by 100.Goodness of fit was assessed for the fitted phenology model using R 2 , standard deviation error in calculation (SDEC), and a visual assessment of residual plots [40,41].We assessed predictive power by separately fitting two validation data sets to the trained model and evaluating R 2 , standard deviation error in prediction (SDEP), and residual plots for the fitted data.We calculated SDEC for the training set and SDEP for the validation sets by dividing the residual sum of squares (RSS) by the number of samples (n) and taking the square root.

Phase II: Mapping of Lateral Shoot Phenology to Model
To acquire lateral shoot data to incorporate into our phenology model, we tracked lateral shoot and leader growth from May 2012-August 2013 at a subset of four sites randomly selected from among our phenology study sites (Figure 1).Each week, we measured the length of the terminal leader and the length of one randomly selected, non-shaded, lateral shoot from the upper three whorls, south side of each of the 25 trees previously selected for inclusion in our phenology study.We began measurements at vegetative budbreak and continued until average leader and lateral shoot length was unchanged from the previous week.
To incorporate lateral shoot phenology into our 2011 phenology model, we first fit the beta sigmoid function to the 2013 lateral shoot data, and then temporally validated that 2013 lateral shoot model using the 2012 lateral shoot data.From the trained and validated model, we derived (1) the time of maximum growth rate (T mLat ), after which reproductive bud initiation occurs; and (2) the time of lateral growth cessation (T eLat ), marking the end of reproductive bud anatomical differentiation.We then mapped T mLat and T eLat to the 2011 phenology model, which allows us to visualize the period of reproductive bud development in terms of leader growth.We followed the same methods to fit the lateral shoot model as we used to fit the 2011 phenology model, including normalizing lateral growth data as percent maximum growth prior to fitting the model.

Phase I: Development of Phenology Model
Temperatures and precipitation were near long-term averages in 2011 [42], when data was collected to train our phenology model of leader growth.In contrast, 2012 was a year of record-breaking temperatures throughout most of the contiguous United States.Average temperatures for the six months from February through July were the highest on record for Michigan and most states east of the Rocky Mountains [42].The extreme temperatures were accompanied by drought at many of our sites, which allowed us to test our phenology model under extreme conditions.The BSF model fit the leader growth data well, and provided a high predictive power both spatially and temporally that was robust to high temperatures and drought (Table 2).SDEC or SDEP were below 0.04 (indicating fitting error ± 4.0%) and R 2 was 0.99 for both the trained model and the model fitted with the 2011 validation data set.Many sites experienced record warm weather in March 2012, resulting in a rapid accumulation of heat sums.Nevertheless, the 2011 phenology model predicted the data well, with SDEP = 0.06 and R 2 = 0.98 for the fitted 2012 validation data set.We found no apparent bias in residual plots from fitted training or test data sets.Model-fitting using the BSF results in parameters that are immediately biologically informative [36].Based on the fitted model, the period of active leader growth ranged from 266 to 1315 • Cd, with the growth rate peaking at 789 • Cd (Table 2; Figure 2a).

Phase II: Mapping of Lateral Shoot Phenology to Model
Consistent with our field observations, we noted differences between model predictions of lateral shoot and terminal leader phenology (Figure 2).Budbreak was slightly earlier (45 • Cd) for lateral shoots, which reached maximum growth rates 94 • Cd ahead of the leader.Growth rate was symmetrical for leader elongation, but the growth rate of lateral shoots dropped rapidly after reaching its maximum.Lateral elongation was complete several days (285 • Cd) before leader elongation.
We obtained a good fit and good predictive power using the BSF to model lateral shoot elongation, with an R 2 of 0.97 and SDEC or SDEP of 0.06 for both the fitted 2013 lateral shoot model and for the 2012 validation data set fit to the 2013 model (Figure 2b; Table 2).We identified key points in lateral shoot phenology to map to our 2011 phenology model (Figure 3): budbreak at 221 • Cd, strobilus initiation shortly after 695 • Cd (T mLat , the time of maximum growth rate), and completion of strobilus differentiation concurrent with growth cessation at 1030 • Cd.By plugging these parameters into our 2011 phenology model, we were able to express shoot phenology in terms of terminal leader phenology.Strobilus initiation began around the time that terminal leader elongation was 35% complete, with differentiation complete once leader growth attained 80% of the maximum length.

Discussion
Using thermal time as the only factor, we developed a phenology model that is able to predict terminal leader and lateral shoot phenology in Fraser fir.The 2011 phenology model was developed and validated using data from nine sites that varied in climate and in site conditions.The predictive power of the model was high even during the extreme heat and drought conditions at the study locations in 2012.This suggests that the model is quite robust, at least within the Upper Midwest, and will predict Fraser fir shoot phenology with a high degree of accuracy, even under the increased temperatures and decreased summer precipitation projected for the region by climate change models [43].
Because reproductive phenology is tied to lateral shoot phenology in Abies [14], our model is able to predict the window during which strobilus initiation and differentiation occur.This window may be used by Christmas tree producers, researchers, and seed orchard managers to manipulate reproductive development through the use of treatments that modify hormonal (e.g., PGRs) or environmental signaling (e.g., temperature and water availability) pathways.These treatments must be applied prior to the period of strobilus initiation (around 695 • Cd, T mLat in our phenology model) and continued until differentiation is complete (1030 • Cd, T eLat ).
The BSF produced a good fit to the shoot growth data and predicted shoot elongation well.However, there are a few limitations to our model.First, it cannot be used to predict shoot growth on individual trees, because it was trained using measurements averaged for each site.Second, it cannot be used to explore differences in vigor among sites, because relative growth was used to fit the model.Third, our model could be improved by the addition of an error term to account for differences among sites.Finally, in every case that we examined, the BSF slightly under-predicted the final shoot length as 99% of the observed maximum length, which may be due in part to underestimating T e .It is unclear whether this error is an artifact of the modeling process (due to the iterative nature of non-linear modeling), or a limitation of the BSF (see [37]).This problem of fit will not affect the utility of our model, but users of any model based on the BSF should be aware of this limitation.In the case of our phenology model, lateral shoot elongation-and therefore strobilus differentiation-may continue slightly past the time predicted by the model, T eLat .
To make our model accessible to Christmas tree growers in Michigan, we have incorporated our model into the Michigan State University Enviroweather platform.Enviroweather provides weather-based IPM and production management tools for agriculture and natural resources [44].In the Fraser fir phenology module within Enviroweather, users can track current shoot emergence based on actual growing degree days and can estimate future shoot growth based on projected GDD accumulation.In this paper, we have provided one example of how our model may be used to guide the management of cone production.Similarly, our model may be adapted to drive other management decisions that are dependent on phenology, such as timing of shearing, control of insect pests and weeds, and frost protection.

Conclusions
We developed a robust model to predict the timing of shoot phenology in Fraser fir.Our model accurately predicted shoot phenology during the unusually hot and dry summer of 2012, indicating that it may serve as a useful tool for projecting shoot phenology, even in the face of climate change.This suggests a role for phenology models in predicting the effects of climate change on tree phenology and growth.The 2011 phenology model predicts the window of strobilus initiation and differentiation, as well as the beginning and end of terminal leader and lateral shoot growth.This information may be used to time practices that reduce cone formation in Christmas tree plantations, or enhance cone production in seed orchards.Our phenology model may provide additional guidance to Christmas tree growers, such as indicating the window for the application of pesticides and predicting when trees will be ready for shearing.The capabilities of this model may be easily extended through the addition of other important phenological events, such as those used to determine variable herbicide application rates throughout the growing season.

Figure 1 .
Figure 1.Locations of Christmas tree farms that cooperated in the 2011-2013 phenology study in Michigan.Coordinates for each site are presented in Table1.

Figure 2 .
Figure 2. Fitted and validated (a) 2011 phenology model of leader growth and (b) 2013 phenology model of lateral shoot growth.Models are fitted using a beta sigmoid function (BSF).The time of maximum growth rate is indicated by T m for terminal leader growth or T mLat for lateral shoot growth.Each symbol indicating training and validation data represents mean growth of 25 trees at one of the measurement sites.Note that predicted shoot growth is the same in every year because the models are based on GDD and not calendar or Julian date.

Figure 3 .
Figure 3. Mapping of Fraser fir lateral shoot phenology to 2011 phenology model of leader growth.T mLat and T eLat bound the period of strobilus bud initiation and differentiation, and indicate the time of maximum growth rate and the end of lateral shoot elongation, respectively.

Table 1 .
Site characteristics of study locations in Michigan, 2011.
[35]erage distance from study plot to weather stations used to obtain growing degree data; superscript indicates number of stations used, if >1; b Web soil survey[34]; c Average annual (June) daily temperature, 1981-2010 U.S. Climate Normals[35]; d Total annual (June) precipitation, 1981-2010 U.S. Climate Normals[35]; † Data were used to train 2011 phenology model.Other sites were used for validation.

Table 2 .
Fitted parameters and measures of fit and predictive power for phenology models of terminal leader growth and lateral growth in Fraser fir in Michigan Christmas tree plantations.
a Standard deviation error in calculation or prediction.Low values indicate good fit to the data; b C m is relative maximum growth rate (% max.growth rate/ • Cd).T b , T e , and T m are in • Cd and indicate the beginning of growth, the end of growth, and the time of maximum rate of growth, respectively.