Willow phenological modelling at different altitudes in central Italy

In order to estimate the impact of climate change on the phenological parameters and to compare them with the historical record, a decision support system (DSS) has been applied employing a Phenological Modelling Platform. Biological observations of two willow species (Salix acutifolia and smithiana Willd) in 3 gardens at different altitudes located in Central Italy were utilized to identify suitable phenological models related to four main vegetative phase timings (BBCH11, BBCH91, BBCH 94, BBCH95), and male full flowering (BBCH 65) clearly identifiable in these species. The present investigation identifies the best phenological models for the main phenophases allowing their practical application as real-time monitoring and plant development prediction tools. Sigmoid model revealed high performances in simulating spring vegetative phases, BBCH11 (First leaves unfolded), and BBCH91 (Shoot and foliage growth completed). Salix acutifolia Willd. development appeared to be more related to temperature amount interpreted by phenological models in comparison to Salix smithiana Willd. above all during spring (BBCH11 and 91), probably due to a different grade of phenotypic plasticity between the 2 considered species.


Introduction
Plant phenological development is mainly influenced by meteorological and climatic factors (Keller 2015), and temperature, solar radiation and precipitation may have a great influence on vegetative and reproductive development, yields and fruit characteristics (Kozlowski and Pallardy 1997;Moriondo et al. 2015;Orlandi et al. 2020). Temperature represents the main factor in the phenological development of plants, which allows the creation of efficient interpretative models for a wide range of species in many different countries and regions (Chuine et al. 2013). Various plant species developments are modelled to predict phenophases and to assess the impact of temperature change on plants (García de Cortázar-Atauri et al. 2010;Ibañez et al. 2010;Pearson 2019;Rojo et al. 2020).
In these studies, the availability of standardized data is essential to assure the replicability of the same monitoring. The "International Phenological Gardens" network through protocols adoption furnish detailed professional observation-guide (http://ipg.hu-berlin.de) and the methodology for following the phenological observations is the same in all the gardens and the references toward the international standardized BBCH method is really recommended. Specific sections of the Phenological Gardens are dedicated to clone plants derived by vegetative propagation through stem cutting of trees and shrubs common to all the network gardens, in order to record only the environmental influence on their development. The first arboreal species utilized in Italian gardens were willows (Salix acutifolia Willd., Salix smithiana Willd. and Salix viminalis L.) which came from central Europe. There is a progressive increase in willow diversity from south to north with the median number of taxa per stand in southern Europe being three, and six in northern Europe (Cronk et al. 2015).
Several models have been used during the last 60 years to simulate phenological development stages; for flowering and fruit setting, the conventional linear growing degree-days model (GDD) was initially proposed by Amerine and Winkler (1944). This simple temperature accumulation model uses daily mean temperatures, with a fixed base temperature and onsets temperature accumulation at fixed dates. Some simulations based on the Winkler model and the Wang and Engel curvilinear model (García de Cortázar-Atauri et al. 2010) were formulated and tested to interpret the flowering phenomenon.
In particular, willows show a strong thermal control over flower buds very sensitive to the first heat amounts during the winter ending period, so GDD accumulations can be very useful to predict reproductive phenology, and the annual flowering cycle reliance on temperature amounts demonstrates the genetic control over phenological response (Mosseler and Papadopol 1989).
Phenological models may also be paired to climate change scenarios so as to project their potential impacts on phenological timings (Ziello et al. 2012;Fu et al. 2015). Therefore, it is critical to understand to what extent temperature may influence the timings of the reproductive and vegetative cycles, as well as identifying plant varieties that better adapt to the climatic changes. The purpose of this study was to determine and analyse the phenological growth stages of two willow species (Salix acutifolia and smithiana Willd.) in comparison to the climate characteristics recorded in three Phenological Gardens of central Italy included in the International Phenological Gardens network (IPG). The phenological records from 2005 to 2018 were used to identify suitable phenological models related to four mains vegetative phenological timings (in spring and fall) and flowering of the two willow species grown in the phenological gardens.

Study area and datasets
Phenological observations from three phenological gardens in central Italy were used for model calibration by monitoring indicator plants left to grow in a Inatural way for as long as possible (Orlandi et al. 2014) (Fig. 1). The first phenological garden, located 15 km from Perugia (Umbria Region, central Italy) in an area characterized by Mediterranean climate, is one of the oldest Italian ones and includes local, national and some indicator species common to all IPGs (Schnelle and Volkert 1964) planted since 1994. The garden has the following geographical coordinates: latitude, 43°00′ 40″ N; longitude, 12°14′ 52″ E, and its total area is of about 4000 m 2 . The second phenological garden considered is that located near Rieti (Lazio Region, central Italy, lat The 2 willow species present in all the phenological gardens and considered by the study were as follows: -Salix acutifolia Willd. Family: Salicaceae; common name: violet willow, sharp-leaf willow; a deciduous shrub or a small tree up to 10 m high flowering during February-April, before the bud burst, and fructifying in May-June. According to what described by Rahmonov (2016) Salix acutifolia seems to be a species that prefers sandy areas typical of humid environments and appears to be a primary colonizing species in its distribution range. -Salix smithiana Willd. Family: Salicaceae; common name: Smith's willow; a deciduous shrub or a small tree up to 9 m high flowering during February-April, before the bud burst, and fructifying in May-June.
The Salix smithiana species is a hybrid species that derives from the crossing of different species of willow such as Salix caprea × viminalis, Salix smithiana for this reason can also be found in our Mediterranean habitats as seen in GBIF Secretariat (https://www.gbif.org).
The plant species were derived several years ago from mother plants received from the German Weather Service as European coordinator for the phenological gardens activities. For model calibration (parameter fitting), temperature records from the meteorological stations installed near the different gardens and over the period of 14 years (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) were considered. The daily meteorological data for the Perugia garden were provided by the station of the Italian Agrometeorological Network (datalogger DA 9000) located about 50 m from the phenological garden. The data for the two Phenological gardens of Rieti and Terminillo Mt. were supplied directly by the meteorological stations located nearby the gardens themselves and managed by the "Apennines Centre of Terminillo Mt." of the University of Perugia (http://www.cat.unipg.it). Furthermore, in order to enhance the statistical robustness of the selected models for a given willow species, the phenological timings for the same species for all three sites were considered for model calibration.

Phenological phases and models
Every biological phase was interpreted by a phenological point of view through observations realized directly in the phenological gardens utilizing the international conventional "BBCH" key to obtain comparable values as reported in literature (Chmielewski and Roetzer 2001;Meier 2001;Saska and Kuzovkina 2010). In the present study, the following phenological phases were considered for the vegetative cycle of both Salix species: BBCH11, first leaves unfolded; BBCH91, shoot growth completed; foliage still fully green; BBCH 94, 50% of leaves discoloured; BBCH95, leaves have mostly fallen (50% of leaves have fallen). Also, BBCH 65 (male full flowering) was reported to highlight reproductive developments during the study period. In each garden, observations were conducted on three individuals of each Salix species, to limit random variability, and the mean date of the onset of each phenophase was calculated as an average of the three plants. Moreover, two principal Leaf presence Periods (LP) were calculated as the differences between the BBCH95-11 dates and BBCH95-91 dates during the study years in the 3 areas. In particular, LP represented the existence of efficient leaves by a physiological point of view considering that phase BBCH 95 is referred to 50% of leaves fallen with the other in a very advanced senescence phase.
An evaluation of the efficacy of the interpretation by the statistical models to follow the biological trends during the yearly LP development (with a step of 5% of the total LP) was realized considering both Salix species in the 3 gardens during the whole study period. The aim of this evaluation was to identify averagely the potential presence of a better simulated LP stage during the entire vegetative season.
Various models were used to simulate thermal accumulation during the considered phenological phases, and the freely adjusted best-fit parameters for different phenological forcing models (GDD, Richardson, Sigmoid and Wang) were selected considering that their performances can be significantly more efficient than those of the default parameters (Caffarra et al. 2011). The best performing models for each phenological phase were selected.

Phenological trend analysis
A trend analysis of the two willows phenological series was realized using nonparametric Mann-Kendall tests, for monotonic positive or negative trends. Positive Z values reveal the presence of a delay trend in the phenological data, while negative Z values indicate an advance trend for the same data. The estimation of a true slope of an existing trend (as change per year) was realized through the Sen nonparametric method. The significance levels were represented by the following symbols: *, trend at p < 0.05; **, trend at p < 0.01; and ***, trend at p < 0.001. The Mann-Kendall trend analysis was realized with the Excel template application MAKESENS version 1.0 (Salmi et al. 2002).

Modelling tools and performance verification
In this study, the Phenological Modelling Platform, PMP (Garcia de Cortazar-Atauri et al. 2010;Chuine et al. 2013;Costa et al. 2019), was used to test and calibrate different phenological models applied to both Salix species. PMP was considered and employed as a useful interactive software-based system (decision support systems, DSS) for decision-making with the use of large volumes of environmental information. The platform allowed to apply and test different phenological models using climatic and phenological data for calibrating the same models and estimating the best-fit model parameters in the specific investigation area. PMP evaluates best-fit parameters on the base of the simulated annealing algorithm of Metropolis et al. (1953) through an iterative optimization procedure.
For each specific location, the environmental and biological data used were daily mean (Td), maximum (Tmax) and minimum (Tmin) temperatures, latitude and observed phenological timings for leaves growth and their senescence and flowering (in days of the year, DOY). The model performance is expressed by the root-mean-square error (RMSE), Nash-Sutcliffe coefficient of efficiency (EF) and then the Akaike Information Criterion (AIC). Particularly, Nash-Sutcliffe efficiency is a normalized statistic that determines the relative magnitude of the residual variance ("noise") compared with the measured data variance ("information") showing the effectiveness of the simulated data reconstruction (Nash and Sutcliffe 1970;Moriasi et al. 2007;Costa et al. 2019). EF varies between -∞ and + 1, + 1 corresponds to an ideal fit, 0 corresponds to the performance of the null model (average), while a negative value to a worse prediction than the null model.

Phenological model performance
The best models for each corresponding vegetative and reproductive phenophases are shown in Table 1 considering at the same time both the Salix species data and the 3 garden areas from 2005 to 2018. EF values evidenced the sigmoid function as the more efficient in the interpretation of all the phenomena in comparison to other models (Caffarra et al. 2011). The T 0 parameters (starting dates of the model calculation for each phase) were always converted to a photoperiod evaluation allowing starting dates to vary with the latitude of the sites. PMP calculated a photoperiod using the latitude value of each site during the calculation period, considering whether T 0 was situated during the increasing day length period (DLP) for spring phenological phases (BBCH 11,91,65) or during the decreasing DLP (BBCH 95) (Caffarra et al. 2011).
Spring vegetative phases (BBCH11 and BBCH91) were those with the best EF values, while the senescence phase (BBCH95) and the reproductive phase (BBCH65) showed lowest values. In particular, full flowering (BBCH65) showed very low values with all the utilized models highlighting how difficult the interpretation of the reproductive process is with the only use of temperature values.
In Table 2, the interpretation of the principal phenophases was realized separately for Salix acutifolia and Salix smithiana utilizing only sigmoid function. Moreover, at first the modelling was realized using only the temperature variable and then considering the Day Length variable; in this case, the Bidabe function was utilized for the second variable. Salix acutifolia showed better results than Salix smithiana for all the phenophases with EF values of 0.79 for BBCH11, 0.75 for BBCH91 and 0.55 for BBCH 95. The full flowering phase did not show sufficient results, as well in neither of the 2 species.
The addition of Day Length in the model functions caused only a small increase in EF values apart from the autumnal senescence phase (BBCH95) which moved from 0.55 to 0.66 for Salix acutifolia and from 0.56 to 0.61 for S.smithiana. Moreover, the model performances are generally higher for acutifolia than for smithiana, which may suggest that the latter is less sensitive to thermal forcing.
In figure 1 the leaf periods (LP) box plots shown have been calculated year by year, as the difference between BBCH95 and BBCH11 dates. The box plots refer to the singular species mean values in the 3 gardens (inter-species phenological effect). The LP median values were around 200 days, considering all the species, while the range of distribution of these data was of 20 days, in average (190-210 days). Moreover, intersite (Perugia, Rieti, Pian di Rosce) phenological effects are shown considering both of the species.
In addition, different LP percentages (with a step of 5% of the total LP) were calculated year by year considering both Salix species in the 3 gardens during the whole study period. The yearly dates (DOY) related to the different LP stages were utilized in the simulation by the modelling platform.
The sigmoid function of the utilized model expressed EF values of the vegetative phases (BBCH11-91-94-95) and of progressive LP percentages that were represented in a unique chart (Fig. 2) consequently during the year.        Fig. 3, as well.
As far as it concerns BBCH11 evidence derived by the Platform, the sigmoid model provided values that were highly conditioned by the strong negative trend manifested during the 14-year study period. Figure 4 a shows an in-depth analysis of the BBCH11 dates during the study period for S. acutifolia, which allowed the detection of a marked discontinuity in the BBCH11 time series, identifying 2 sub-periods, BBCH11A, from 2005 to 2010 and BBCH11B from 2011 to 2018, characterized by similar regression lines (Fig. 4b).

Discussion
Phenological models may produce valuable information for the fields of agriculture and forestry, suggesting better practices and the implementation of appropriate measures to optimize yields (Cola et al. 2014;Rojo and Pérez-Badia 2015). The timing of yearly leaves presence from first leaves unfolded, to shoot growth and fallen leaves, influences the photosynthesis, affecting biomass production, but can also expose plants to the risk of spring late frosts (Polgar and Primack 2011). The leaf development phenomenon is very critical for modelling present and future forest and tree-crop dynamics, since it can determine the availability of food and shelter for many species, creating cascading effects through multiple trophic levels (Wesolowski and Rowinski 2006;Both et al. 2009;Coyle et al. 2010). The present study allowed the identification of the best phenological models in willow species for which Sigmoid model performed efficiently in simulating spring vegetative phases, BBCH11 and 91 respectively. However, for flowering the performances were not satisfactory probably due to the fact that in deciduous arboreal species with precocious flowering dates at the end of winter-early spring, the reproductive development is particularly influenced by the chilling amounts (flowering induction phenomena) in addition to the successive forcing temperature (after break dormancy). This winter chilling requirement clearly represents an adaptation of the plant in order to prevent the beginning of reproductive development may the frosts during February-March in the Mediterranean area become very dangerous (Larcher 2003).
As regard the calculation of the leaf period (LP) and its subdivision in progressive percentages during the vegetative development, the phenological platform averagely designated the month of June, when full leaf development (BBCH91) was recorded during the study period, as the period that was better simulated by the models. This is probably because this is the phase of the year when the forcing temperatures express their influence mainly on vegetative structures, while precedent phenological phases could be influenced also by winter chilling effects. Successive phenophases such as the senescence ones during September-October could be influenced by summer water availability stress, therefore have a reduced efficiency of model interpretation.
The modelling results showed as the addition of Day Length variable in the interpretative models was significant above all during the senescence phase when the light hours decrease rapidly confirming as photoperiod play an important role in regulating the leaf development of temperate woody plants (Ghelardini et al. 2014). Decreasing day-length induces growth cessation and bud set in many perennial plants (Lagercrantz 2009) and, along with temperature, it may force leaves senescence and shedding (Kozlowski and Pallardy 1997;Tanino et al. 2010;Heide 2011). Several studies indicate that temperature and stress factors may variably interact with photoperiod in controlling the timings of phenological events in woody species, including willow and poplar (Barros and Neill 1987;Molmann et al. 2005;Kalcsits et al. 2009).
As regard the trend analysis results, the BBCH11 time series discontinuity, with the definition of two sub-periods, was due to several positive winter temperature anomalies recorded from 2011 to 2018. In this last 8-year period, only winter 2012 showed a negative temperature anomaly (limited to a polar jet stream phenomenon during 2 weeks of February). All other winters during this period appeared to be averagely milder in the  (Merilä and Hendry 2014;Urban et al. 2016). In temperate plants, the leaf emission phenomenon may also be highly adaptable. In this view, while species with low plasticity and adaptability may suffer under warmer conditions, other species will likely thrive. However, the same phenotypic variation effect is actually discussed mainly in presence of low specific genetic variation, considering that, in the short term, the only way plants can respond to temperatures rising is by adjusting to the new environmental conditions through plasticity, while, in the long term, they achieve this through adaptive evolution, by exploiting both the presence of genetic variation and the seasonal adaptation (Oostra et al. 2018).
Salix acutifolia showed phenological development stages more efficiently interpreted through simulation by the utilized statistical models in comparison to Salix smithiana, especially during spring (BBCH11 and 91). In this sense, S. acutifolia showed highly temperaturemediated phenological traits probably due to a specific phenotypic plasticity that will ensure the survival of the species populations under future climate change projections, even if generalizations are very difficult considering that some species are limited by temperature and photoperiod (Savolainen et al. 2011). In the future  climate scenario, the geographical area of distribution of Salicaceae plants is generally likely to shift to high latitudes. Salicaceae is a group of plants with a large ecological amplitude, thanks to which they can adapt to many different ecological environments. However, if global warming becomes more intense, the suitable range of Salicaceae faces the risk of reduction. Species of Salicaceae are, therefore, clearly vulnerable to climate change effects, with contractions in their range being mainly concentrated at low latitudes (Bertrand et al. 2011;Garcia et al. 2013). Lastly, the phenological derived models will be used to assess the impacts of climate change on the critical development phases in different species, running simulations with the use of climatic data from different forcing scenarios and climate model experiments. Recent simulations of the potential changes in phenological timings may also deliver important information for medium-to-long term planning, suggesting future anticipations of the spring phenophase timings up to about 30 days until the middle of the twenty-first century (Ibanez et al. 2010;Pearson 2019). The growing season length study and forecasting can have important impacts on ecosystem processes, including the uptake of CO 2 , tree growth, microclimate and water movement (White et al. 1999;Morisette et al. 2009). Tree foliage investigations may suggest the ecological consequences of such growing season increase in terms of carbon sequestrations balance due to a contemporary variation of two biological phenomena, photosynthesis and respiration inside forests (Milyukova et al. 2002;Dunn et al. 2007;Piao et al. 2008;Richardson et al. 2010).
Funding Open access funding provided by Università degli Studi di Perugia within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.