Modeling phenological responses of Inner Mongolia grassland species to regional climate change

Plant phenology is an important indicator of ecosystem dynamics and services. However, little is understood of its responses to climate change, particularly in ecologically sensitive regions such as arid and semi-arid grasslands. In this study, we analyzed a long-term climate and plant phenology dataset of thirteen grassland species in the Inner Mongolia of China, collected during 1981–2011 time period, to understand temporal patterns of plant phenology and then developed a simple chilling-adjusted physiological model to simulate phenological responses of each plant species to climate change. The results of regional climate analysis suggested that the minimum temperature was increasing at a greater rate than mean and maximum temperatures in the region and the climate variability had significant impacts on vegetation phenology. Chilling from an early stage in spring in general slowed down the phenological development in most plant species, although there were some inconsistencies among sites and years. Specifically, we found lower precipitation during green-up resulted in delayed flowering, which may attribute to plant self-adjustment strategy to respond changes in climate. These climate dependent phenologies were characterized by a simple physiological model. Scenario analysis suggested that by 2071–2100 significant shifts in plant phenology are expected in Inner Mongolia, including as much as 6–11 days earlier in green-up time and 8–11 days shorter in growing season due to earlier senescence.


Introduction
Plant phenology, the timing of recurrent biological events of plants, is determined by biotic and abiotic factors, and by the interrelation among phases of the same or different species (Lieth 1974). The timing of phenological events is a major determinant of plant productivity and species distribution (Vitasse et al 2011). It is related to seasonal and inter-annual dynamics of vegetation and is highly sensitive to climate variability , Han et al 2015, Keenan and Richardson 2015. Therefore, phenological characteristics can be used as an integrated indicator of climate change.
Changes in phenological phases can potentially impact terrestrial ecosystem functions and services such as carbon, water and nutrient cycles, and productivities (Schwartz 1998, White et al 1999, Menzel 2000, Piao et al 2007. Changes in these critical ecosystem services can have significant implications to the livelihoods of human society, particularly rural families such as ranchers and farmers as their lives rely on these services for sustainable development. Plant phenology of grasslands is also a major factor affecting plant function such as light and water use efficiency, seed production, and the total biomass production (Vitasse et al 2011). Climate change can lead to a failure in flowering or fruiting  and result in a low seed production and even a change in the structure of plant communities (Lesica and Kittelson 2010). However, knowledge on how climate change affects phenology of different grassland species is very limited.
In addition to physiology, plant phenology is also dependent on the external environmental temperature (Thompson and Clark 2006). The timing of a given phenological stage is often determined by a threshold of accumulated temperature, which is usually measured by growing degree-days or GDD (McMaster and Wilhelm 1997, Bonhomme 2000, Cook et al 2010. Plant development can be characterized either by a linear function of temperature when the temperature is above a baseline value or by a nonlinear function when it is below the baseline to account for the base, optimal and maximum temperatures for growth. There is also a negative impact term accounting for extreme high temperature (i.e. supra-optimal temperature) that may progressively deactivate certain plant enzymes, and consequently slow down growth (Atkinson andPorter 1996, Bonhomme 2000).
Previous studies revealed that a period of low temperature (chilling effect) during an early development stage could hasten flowering of some herbaceous (e.g., Cook et al 2012 and some woody species (e.g., Luedeling et al 2009, Richardson et al 2013, Clark et al 2014. Sometimes it can also promote flowering of some perennial grasses in temperate and cold environments (Heide 1994, Atkinson and Porter 1996, Shirazi 2003. As global warming trend continues, especially in winter times, the chilling period may be shortened or even disappear in some geographic regions to affect early vegetation growth (Cook et Guo et al 2015). Shortening or disappearance of chilling period, due to climate warming, has significant implications for grassland species that require chilling effect to trigger growth (Donald 1984), such as Kobresia pygmaea as shown in a warming trend experiment by Dorji et al (2013). To date, however, little is known about chilling effect on phenological events of individual grassland species in Mongolian plateau.
It is imperative to understand phenological change at individual species level in grasslands in response to seasonal and daily temperature variations, particularly to chilling effect, which has significant implications in species distribution, biomass production, and water and carbon cycles, in anticipation of future climate change. The objectives of this study are to: (1) characterize chilling effect on plant phenology at species level using a long-term phenological dataset from Inner Mongolia, China; (2) develop a physiological model to account for and subsequently simulate chilling effects on species phenology; and (3) predict future climate change impacts on species distribution in Inner Mongolian grasslands.

Study sites
Phenological and climate data were collected at eight (8) experimental sites in Inner Mongolia Autonomous Region of China, located between 97°10′E to 126°09′E and 37°24′N to 53°20′N (figure 1). The vegetation types at these sites include meadow steppe (found in Ergun, Ewenkizu Qi or Ewenki, and Bayaertuhushuo or Bayar), typical steppe (found in Xilinhot or Xilinhot, Xianghuang Qi or Xianghuang and Qahar Youyi Houqi or Qahar) and desert steppe (found in Wushenzhao or Wushenzhao and Alxa League or Alxa). The climate in Inner Mongolia can be characterized by cold and dry winter and warm and rainy summer, with higher rain-use efficiency than most other areas in the semi-arid and arid regions of China (Zhang et al 2011). The annual rainfall decreases gradually from the East (∼400 mm) to the West (∼200 mm) while the corresponding annual mean temperature ranges from −1.6°C to 9.2°C.

Phenological and climate data
The phenological data of thirteen dominant grassland species were collected at the eight experimental sites, including the timing of the main phenophases, such as green-up, heading, flowering, fruiting and senescence from 1981 to 2011 by the China Meteorological Administration and the Inner Mongolia Meteorological Bureau of China. These experimental sites were located in natural pastures with an area of ∼10 000 m 2 at each site where thirteen dominant grassland species (table 1) were found to include Poaceae (6 species), Rosaceae (1 species), Asteraceae (2 species), Zygophyllaceae (2 species), Fabaceae (1 species) and Tamaricaceae (1 species). A permanent experimental plot (10 m 2 ) at each site was set up to collect phenology data for each species found within the plot. Ten individual plants for each species were examined in the afternoon every two days to record phenological information.
The date of each phenological event was recorded when 50% of observed plants reached to a phenological event. For species with branches, the observations were made on the main stems only.
In addition to the phenological data, daily minimum, mean and maximum air temperatures along with daily precipitation data from 1981 through 2011 were obtained from the China Meteorological Administration recorded at nearby weather stations (supplementary table 1).
For future scenario analysis, climate data from 2041 to 2100 timeframe were taken from the output of PRECIS (Providing Regional Climates for Impacts Studies), a regional climate modeling system based on the third generation of the Hadley Centre's regional climate model (HadCM3) (Jones et al 2004). The IPCC SRES A2 scenario, which represents an environmentally pessimistic development pathway in China under medium-high emissions, is used in this study. The model predictions from 1961 to 1990 were in a good agreement with the climate record in Inner Mongolia (Xu et al 2006) and were used in estimating the grassland net primary productivity . The predicted future phenologies of the species were averaged from 2041 to 2070 to present the period of the 2050s (s indicates a decade), and from 2071 to 2100 to represent the period of the 2080s. The future phenological characteristics were compared to the historical observation data from 1981 to 2011 to derive standard deviations of phenology in the future. In this study, the observed phenological events and climate data from 1981 to 2000 were used for model development, while the data from 2001 to 2011 were used for model verification. Since the observation data of Caragana stenophylla at Alxa was only available from 1983 to 2002, the data prior to 1995 for this species were used for model development and those acquired after 1995 were used for model validation.

Description of phenological models
In this study, we first analyzed phenological data from 1981 to 2000 using an existing thermal time (TT) model and an existing physiological time (PT) model. The PT model then was modified to include chilling effect (PTc) to assess potential impact of climate warming trend on plant phenologies. The new PT C model was then validated using data from 2001 to 2011 and subsequently used to simulate phenologies under future climate scenarios.
The TT model is based on the concept of TT (°C d), similar to models based on growing degree days (Monteith 1977). It assumes a linear relationship between development rate and temperature above a base line T b (Bonhomme 2000, Thompson andClark 2006, Zhang et al 2008): where TT is set to zero when the temperature T is below T b . The PT model assumes a 'broken stick' response function (see supplementary figure 1) with the temperature (Diekmann 1996, Campbell and Norman 1998, Soltani et al 2006, Zhang et al 2008. It further assumes a maximum development rate when temperature is above a threshold and it calculates the number of days required for a certain phenological event to occur (Zhang et al 2008. The PT can be either calculated from field observations or modeled as a function of temperature. Previous studies suggested that the function is usually a two-phase linear 'broken stick' or a smooth optimum relationship with temperature (Yin et al 1995).
It is noted that in this study we did not consider a maximum temperature tolerance (e.g. 35°C for most plant species) that could hinder plant development, because daily temperature during the warmest month, July, in Inner Mongolia rarely reaches 23°C. Therefore, we assume that the daily growth is at its maximum rate when the temperature is above the optimum baseline.
In this study, PT is calculated by summing daily values of the thermal effectiveness f T ( ) (ranging from 0.0 to 1.0) starting from the first day of the year (Hanninen 1995, Cao andMoss 1997): T is the daily temperature, T b is the base temperature, and T o is the optimal temperature.
In addition to the accumulation of heat needed for plants to develop leaves or flower in spring, many plants in temperate climate zone require an accumulation of cold days (chilling effect) to break plant dormancy for woody species (Campoy et al 2012, Luedeling et al 2013, Guo et al 2015 and to vernalize floral initiation in herbaceous species (Donald 1984, Cao andMoss 1997). Therefore, the PT model was further modified to account for chilling effect by assuming that lack of chilling would delay the optimal development. The modified model, PT C, accounts for both the effect of chilling completion (CC) and thermal effectiveness to achieve phenological stages.
The chilling effect, c(T), can be modeled by setting the value to be 1.0 when the temperature is between the baseline and optimal temperatures (Cao and Moss 1997) and 0.0 when the temperature is below chilling-base temperature T cb to account for frozen effect. When the temperature is in between T co and T cm , it is a simple linear function. There is no chilling effect when the temperature is above the maximum chilling temperature (T cm ) or below the chilling-base temperature (T cb ). Mathematically chilling effect, c(T), can be expressed as: where, T cb is the base temperature to be determined experimentally, T co is the optimal temperature, and T cm is the maximum temperature for chilling effect. Phenological development is dependent on CC (Cao and Moss 1997), expressed as the fraction of accumulated chilling effect c(T) from 1 October of the previous year to the total chilling requirement (CC m ). The value of CC is set to a maximum value of 1.0 Chilling-adjusted PT (PT C ) can be modeled as the product of daily temperature effect f(T) and CC (Cao and Moss 1997). The PT C is the number of physiological days required under optimal development and CC. If the chilling is not completed, the development rate of species is reduced proportionally. Due to the chilling requirement of floral initiation in grass species, the chilling effect is only used for predicting flowering events. Therefore, the PT c is expressed as: The integration time (t, days) of TT, PT and PT C starts from 1 January of a year to each event for all species. However, the CC is calculated from 1 October of the previous year to include both initiations of leaf and the floral buds (Donald 1984, Chen et al 2014.

Model parameterization
The base temperature (T b ) and optimal temperature (T o ) for the model were determined by fitting the model to historical observations. Doing so, we first tested T b using a TT model (Monteith 1977, Bonhomme 2000, Thompson and Clark 2006, Zhang et al 2008, by setting T b value from 0°C to 10°C to determine the T b with the least square fitting. We then tested T o by using the PT model at a fixed T b determined in the first step with the value of T o set between 10°C and 30°C to determine the T o with least square fitting. To determine the threshold temperatures of chilling, we first fixed the chilling-base temperature for grass species based on previous research , and then determined the values of optimal and maximum chilling temperatures. Further, based on the analysis of experimental data, T b and T o were set to 0°C and 15°C for equation (3), respectively. T , cb T co and T cm were set to −5°C, 0°C, and 5°C for equation (4), respectively. The same thresholds were used in the model for all herbaceous species and shrubs found at our study sites in Inner Mongolia.
The requirements of CC m for different species were categorized into strong (CC m =65), moderate (CC m =60) and weak (CC m =50) sensitivities to chilling temperatures. The CC m was estimated from the phenological data from 1981 to 2000.

Accounting for diurnal temperature variation
The daily maximum and minimum temperatures do not change at the same rate; they are often accompanied by reduced frost occurrence and increased heat stresses . Including the variation of minimum and maximum temperatures rather than only using daily average in the model enables better presentation of phenological characteristics in understanding climate impacts on plant growth.
The daily mean temperature was calculated as the average of the daily maximum and minimum values (Cannell and Smith 1983). The effect of diurnal temperature variation on chilling-adjusted PT (PT C ) was considered as 50% in average, 25% in both maximum and minimum air temperatures (Yang et al 2004, Zhang et al 2008. The actual daily thermal effect f(T) can be distributed among these three temperatures by: here f(T mean ), f(T max ) and f(T min ) are the thermal effect calculated in equation (3), replacing T with T mean , T max and T min , respectively. Calculation of actual c(T) regarding daily temperature variation is similar to the calculation of actual f(T).

Historical phenology and model evaluation
To assess the predictive capability of the model, the root mean square error (RMSE) and the coefficient of variation of the root mean square error (CVRMSE) were used as indicators of model fitting:  2(a)). In other words, minimum temperatures were increasing at a rate that is faster than mean and maximum temperatures. Annual precipitation indicate an insignificant decreasing trend (1.53 mm per decade) (R 2 =0.08, P=0.1) ( figure 2(b)).

Phenological thermal requirements
The thermal requirement for each phenological development stage differed significantly (P<0.01) among species (  (table 3), which indicates the possibility of difference in CC for different species with respect to chilling sensitivity (strong, medium or weak). The chilling time in meadow steppe or typical steppe was less than 50 while in desert steppe, it was slightly more than 60. The variation (SE/mean, %) in chilling-adjusted PT (PT C ) was 2.42% for senescence. The variations of PT C for all developmental stages were similar to that of PT and smaller than that of TT. These results suggest that all species in Inner Mongolia meet the chilling requirement under current climatic condition.

Model validation
Compared with the phenology observations from 2001 to 2011, the PT C model results showed agreement with observations for the timing of green-up, heading, flowering, and fruiting. The overall prediction bias across all species, years, and sites was −2.4 days for time to green-up, 3.4 days for flowering, −1.7 days for fruiting, and −7.6 days for senescence. The overall CVRMSE was 9.0% for green-up, 6.9% for flowering, 6.3% for fruiting, and 6.8% for senescence (table 4). For all phenological events, the chillingadjusted PT (PT C ) models for Stipa and Leymus chinensis had small biases ( figure 3). However, the predictions at the Xianghuang site showed a large discrepancy compared to the observations (figure 3), which we postulate was caused by frequent droughts.

Simulated phenological change under future climate change
Using the A2 climate scenario, we projected phenology for Stipa and Leymus chinensis in two future time periods: 2041-2070 (2050s) and 2071-2100 (2080s). For Stipa, the timing of green-up in the 2050s is 6.4 days earlier than in the 2080s, but not significantly different from historical observations. The length of the growing season of Stipa would be shortened by 9.8 days in the 2050s and by 11.0 days in the 2080s (figure 4(a)). For Leymus chinensis, the green-up time was 5.8 days earlier in the 2050s and 10.9 days earlier in the 2080s than historically observed values. The length of the growing season of Leymus chinensis is predicted to be shortened by 9.0 days in the 2050s and by 7.8 days in the 2080s ( figure 4(b)). The shortening of growing season for two grass species in future climate change comparing to historical observations was due to earlier senescence.

Discussion and conclusions
The plant phenology in Inner Mongolia varied between years and among the study sites due to difference in climate patterns. The species in desert steppe in the Western part of the region required much more TT (°C d) or PT (days) than those in meadow and typical steppe in the central and Eastern part of the region.  The chilling-adjusted physiological phenology model captured fundamental characteristics of species, indicated by its agreement with the observations. The bias for each modeled phenophase ranged from −7.6 to 3.4 days with CVRMSE<9.0%. It should be noted that the current model neglects photoperiod Table 3. Chilling-adjusted physiological time (PT C , days) of 13 species at eight observation sites from 1981 to 2000. The mean values with s.e. (standard error for years) was given. All development stages started from the first day of a year. Chilling effect of previous year was included in the model.

Species
Site  effect (Cao and Moss 1997), which is also important in grass species. The chilling-adjusted phenological model did not improve the prediction much in comparison with physiological model (see supplementary  table 2 and figure 2). However, it is more adaptable to assess climate change, particularly when the increase in minimal temperature might reduce the possibility of CC. Compared with the historical records, phenological development from green-up to senescence of two major grass species is projected to be shortened by 9 to 10 days in the 2050s and 8 to 11 days in the 2080s under A2 future climate scenario.
The plant thermal requirement is an essential indication to estimate the length of different phenological development phases (McMaster and Wilhelm 1997, Bonhomme 2000, Trudgill et al 2005. The linear relationship between development rate and temperature is commonly used, and the nonlinear relationship involving optimal and maximal temperatures is   increasingly applied to improve the estimation accuracy as it accounts better for climate extremes (Wang 1960, McMaster and Wilhelm 1997, Agonga 2006. In this study, we used 0°C as a base temperature and 15°C as an optimal temperature for all species based on historical observations. However, controlled experiments are needed to obtain more accurate estimates of cardinal temperatures for different species. Extremely high temperatures would delay plant development rate by deactivating certain enzymes (Bonhomme 2000) or damaging plant tissue (Alshareef 2010). In this study, the negative effect of high temperature was not considered because this region is mostly below optimal temperature. Unlike many woody species where the flower bud is exposed to low temperatures during autumn and winter, the grass species starts floral stage after the physiological grass growth in spring. However, the chilling effect on leaf buds of grass species might start from the autumn in previous year; the accumulation of chilling is therefore reasonable in this study to be set from 1 October of the previous year. Most perennial grasses in temperate climate have a chilling requirement for flowering, with a primary induction brought about by low temperatures (Heide 1994, Tooke andBattey 2010). As global warming trend continues, especially in winter, plant development in spring is expected to be delayed due to the lack of chilling effect (Yu et al 2010, Guo et al 2015. It is difficult to quantify the collective warming and chilling effects on the timing of phenophase. Further, decoupling of chilling followed by climate warming is also difficult to evaluate, with little evidence of chilling effect being found (Clark et al 2014).
There are many practical questions that remain unanswered including if species adapt to warming conditions by changing their chilling requirements, but they should be addressed in future studies (Schwartz and Hanes 2010). This study only explored chilling effects on grass phenology at the species level. Some herbaceous species, especially those that are sensitive to chilling, such as winter wheat, requires certain accumulation of chilling temperature to initiate the flowering (Donald 1984) and thus, chilling becomes an important factor in the phenological model (Cao and Moss 1997). Chilling effect on grass species is typically ignored in many phenological models because its leaf phenology is likely opportunistic. Grass species cannot be vernalized as imbibed seed (a stage when the seeds absorb water from the soil before germination and emergence occur) and therefore must first produce a certain amount of vegetative growth as a seedling, e.g. reed canary grass (Pbalaris Arundinacea) (Donald 1984). The chilling-adjusted phenological model proposed in this study made it possible to study phenological changes in other ecological regions.
The photoperiod effect is also important to predict grass species phenology (Tooke and Battey 2010). However, this is not included in our current phenological model because our observations were made under an environmental condition where the differences in day length were too small to quantify this effect. Photoperiod is found to compensate for the lack of chilling effect, and might be more influential than temperature in some tree species phenology (Way and Montgomery 2014). However, we do not know if this would be the case for leaf phenology (versus flower phenology) of grass species.
Precipitation is a dominant factor determining productivity of grasslands in arid and semi-arid regions  and also an important regulator of spring phenology (Cong et al 2013). The influence of precipitation on phenological phases is also vegetation type and soil water dependent (Jolly and Running 2004, Piao et al 2006, Liancourt et al 2012. We found the flowering and senescence times of Stipa and Leymus chinensis were significantly sensitive to drought (figure 5), but there was little impact on green-up (figures 5(a) and (b)).
There was also an overall increasing trend for precipitation and physiological development for Leymus chinensis, with a regression of y=0.13x+42.92 (R 2 =0.24, P<0.01) ( figure 5(d)). To further improve the prediction of phenology in grass species, patterns of precipitation need to be taken into account (Cong et al 2013, Chen et al 2014, Shen et al 2014. In the future, it would be beneficial to quantify the drought effect on grass phenology while considering both regional precipitation and soil moisture content (pre-green-up precipitation). Additionally, plant fertilization is also sensitive to temperature stress during the flowering and pollen development stages (Thuzar et al 2010).
Change in phenological events and growing season length, associated with the future climate change, could be crucial for determining appropriate grazing pressure, estimating net primary productivity of grasslands (Piao et al 2007, Caldararu et al 2014, Chen et al 2014, and ecosystem functioning (Edwards and Richardson 2004). Controlled experiments are needed to understand the physiological and ecological mechanisms on the relationship between species-level phenological events and thermal requirement, such as chilling effect.