Bertold Mariën

Manuela Balzarolo

Inge Dox

Sebastien Leys

Lorène J Marchand

Charly Geron

Miguel Portillo-Estrada

Hamada Abdelgawad

Han Asard

Matteo Campioli

Mariën Bertold

Balzarolo Manuela

Dox Inge

Leys Sebastien

Marchand J Lorène

Geron Charly

Portillo-Estrada Miguel

Abdelgawad Hamada

Asard Han

Campioli Matteo

Orcid Id

Bertold Marien
email: bertold.marien@uantwerpen.be

Hamada Abd Elgawad

Detecting the onset of autumn leaf senescence in deciduous forest trees of the temperate zone

Keywords: Phenology, Autumn senescence, Breakpoints, Canopy dynamics, Fagus sylvatica, Betula pendula, Quercus robur, Leaf coloration and fall

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Autumn leaf senescence is a controlled type of programmed cell death that, unlike other stressors inducing cell death, avoids the loss of leaf nutrients [START_REF] Keskitalo | A cellular timetable of autumn senescence[END_REF]. While annual plants use the senescence process to transfer nutrients from their leaves to their maturing seeds, perennials specifically recapture leaf nutrients during autumn to relocate them to their over-wintering organs, as this is essential for their growth potential and foliage redevelopment during the subsequent year [START_REF] Hagen-Thorn | Autumn nutrient resorption and losses in four deciduous forest tree species[END_REF][START_REF] Weih | Genetic and environmental variation in spring and autumn phenology of biomass willows (Salix spp.): effects on shoot growth and nitrogen economy[END_REF][START_REF] Estiarte | Alteration of the phenology of leaf senescence and fall in winter deciduous species by climate change: effects on nutrient proficiency[END_REF]

Accepted Article

. [START_REF] Estiarte | Alteration of the phenology of leaf senescence and fall in winter deciduous species by climate change: effects on nutrient proficiency[END_REF]. For annuals, leaf senescence has been studied in details because of its relevant applications in agriculture and food storage [START_REF] Buchanan-Wollaston | The molecular analysis of leaf senescence--a genomics approach[END_REF][START_REF] Diaz | Characterization of markers to determine the extent and variability of leaf senescence in Arabidopsis. A metabolic profiling approach[END_REF]. Ecophysiological and, in particular, molecular approaches have elucidated the main aspects of this process, not only describing the molecular pathways related to leaf senescence, but also its potential drivers. For instance, cytokinins have been found to play a key role in annuals by blocking or delaying leaf senescence [START_REF] Kudo | Metabolism and long-distance translocation of cytokinins[END_REF][START_REF] Edlund | Contrasting patterns of cytokinins between years in senescing aspen leaves[END_REF].

For deciduous trees, the situation is different and we are far from understanding the molecular mechanisms controlling leaf senescence in this plant functional type. However, research on Populus tremula identified the timetable of the leaf senescence process [START_REF] Andersson | A transcriptional timetable of autumn senescence[END_REF][START_REF] Keskitalo | A cellular timetable of autumn senescence[END_REF][START_REF] Fracheboud | The control of autumn senescence in European aspen[END_REF]. One of the first key stages taking place during leaf senescence is chlorophyll degradation, which leaves undergo to recover valuable elements -especially nitrogen and magnesium-contained in their pigments. During the controlled dismantling of the photosynthetic apparatus, molecules are broken down and conversed into transportable compounds through catabolic enzymes [START_REF] Matile | Biochemistry of Indian summer: physiology of autumnal leaf coloration[END_REF][START_REF] Andersson | A transcriptional timetable of autumn senescence[END_REF][START_REF] Keskitalo | A cellular timetable of autumn senescence[END_REF]. The energy that fuels this relocation phase is first provided by photosynthesis and then, at a later stage, by leaf nonstructural carbohydrates through mitochondrial respiration [START_REF] Keskitalo | A cellular timetable of autumn senescence[END_REF][START_REF] Fracheboud | The control of autumn senescence in European aspen[END_REF]. The chlorophyll degradation does progressively lead to the leaf coloration as other non-green pigments, such as carotenoids, become visible [START_REF] Feild | Why leaves turn red in autumn. The role of anthocyanins in senescing leaves of red-osier dogwood[END_REF][START_REF] Croft | Leaf Pigment Content. Reference Module in Earth Systems and Environmental Sciences[END_REF]. At the end of the leaf senescence process, nutrient transport through the phloem is stopped by the formation of an abscission and separation layer in the petiole, which leads to abscission of the leaf [START_REF] Keskitalo | A cellular timetable of autumn senescence[END_REF][START_REF] Ruttink | A molecular timetable for apical bud formation and dormancy induction in poplar[END_REF][START_REF] Fracheboud | The control of autumn senescence in European aspen[END_REF].

Based on the timetable described above, it is clear that the leaf senescence process comprises different phases which can be characterized by multiple phenological events [START_REF] Gallinat | Autumn, the neglected season in climate change research[END_REF]. One of the key phenological events of the leaf senescence process is the onset of leaf senescence [START_REF] Keskitalo | A cellular timetable of autumn senescence[END_REF][START_REF] Fracheboud | The control of autumn senescence in European aspen[END_REF]. The onset of leaf senescence coincides with the start of chlorophyll degradation and nutrient relocation, and, once in place, it cannot be reversed. The onset of leaf senescence has been defined and investigated by fundamental physiological studies measuring chlorophyll degradation and nutrient relocation but seldom determined using statistics [START_REF] Fracheboud | The control of autumn senescence in European aspen[END_REF][START_REF] Edlund | Contrasting patterns of cytokinins between years in senescing aspen leaves[END_REF]. On the other hand, ecological studies and large scale monitoring programs regularly infer the timing

Accepted Article

. of the leaf senescence process using the later stages of the process (e.g. leaf abscission), which are easy to assess but are physiologically less relevant than the onset of the leaf senescence [START_REF] Gallinat | Autumn, the neglected season in climate change research[END_REF][START_REF] Gill | Changes in autumn senescence in northern hemisphere deciduous trees: a metaanalysis of autumn phenology studies[END_REF].

Assessments of leaf senescence onset are commonly based on three processes taking place during leaf senescence: (1) leaf nutrient relocation, (2) chlorophyll degradation and (3) leaf coloration (Table 1). Measurements of the temporal evolution of the first two processes represent the most direct way to detect autumn leaf senescence. The first method derives onset of leaf senescence from the start of a decline in the content of specific nutrients such as nitrogen (N) in the leaves [START_REF] Keskitalo | A cellular timetable of autumn senescence[END_REF][START_REF] Homolová | Review of optical-based remote sensing for plant trait mapping[END_REF]. The second method derives onset of leaf senescence from the start of a decline in leaf chlorophyll content [START_REF] Keskitalo | A cellular timetable of autumn senescence[END_REF][START_REF] Fracheboud | The control of autumn senescence in European aspen[END_REF]. The main issue associated with these methods is the determination of the exact starting point of the seasonal decline. Statistical methods typically used to solve similar problems require high-frequency time series not easily available from the leaf N and chlorophyll content measurement methods [START_REF] Liu | Using canopy greenness index to identify leaf ecophysiological traits during the foliar senescence in an oak forest[END_REF]. Therefore, a common solution is to take a threshold in the loss of leaf N or chlorophyll (e.g. 20-50% loss compared to the summer maximum) as the onset of leaf senescence instead of the start of the seasonal chlorophyll decline [START_REF] Gunderson | Forest phenology and a warmer climate -growing season extension in relation to climatic provenance[END_REF]. Yet, this is only an approximation. Other issues related to these methods are their laborious character (e.g. dry combustion of ground samples for leaf N, extraction and spectrophotometric analysis for chlorophyll) and potential difficult sampling (e.g. sampling of leaves at the top of the canopy requires a scaffolding tower, tree climbers etc.). The chlorophyll based method can be made simpler by using indirect proxies to measure chlorophyll (e.g. using a chlorophyll content meter) or, when studying forest stands, remote sensing approaches, which can avoid the leaf sampling phase and can make use of several established remote sensing indices of seasonal chlorophyll trends [START_REF] Richardson | An evaluation of noninvasive methods to estimate foliar chlorophyll content[END_REF][START_REF] Dash | The MERIS terrestrial chlorophyll index[END_REF]. On the other hand, no suited direct remote sensing index currently exists for leaf N [START_REF] Homolová | Review of optical-based remote sensing for plant trait mapping[END_REF][START_REF] Hawryło | Estimating defoliation of Scots pine stands using machine learning methods and vegetation indices of Sentinel-2[END_REF]. Following the third approach, the onset of leaf senescence can be derived from the loss of canopy greenness, a variable that considers primarily coloration but also leaf loss of the colored leaves [START_REF] Vitasse | Responses of canopy duration to temperature changes in four temperate tree species: relative contributions of spring and autumn leaf phenology[END_REF]. This approach involves very simple visual observations and is therefore widely used [START_REF] Gill | Changes in autumn senescence in northern hemisphere deciduous trees: a metaanalysis of autumn phenology studies[END_REF]. However, visual observations can be subjective, and the lag introduced (leaf coloration occurs only after a certain degree of chlorophyll degradation and nutrient relocation has taken place) is of uncertain magnitude [START_REF] Christ | Mechanism and Significance of Chlorophyll Breakdown[END_REF]. As the other methods, for forest stands,

Accepted Article

. remote sensing indices related to autumn coloration can be used to standardize the methodology. The plant senescence reflectance index (PSRI; a proxy specifically designed to follow leaf senescence), or coloration proxies such as relative green-or redness, are some examples [START_REF] Chipkin | Effect of propranolol on antinociceptive and withdrawal characteristics of morphine[END_REF]Merzlyak et al., 1999;[START_REF] Sonnentag | Digital repeat photography for phenological research in forest ecosystems[END_REF][START_REF] Filippa | Phenopix: A R package for image-based vegetation phenology[END_REF][START_REF] Richardson | Tracking vegetation phenology across diverse North American biomes using PhenoCam imagery[END_REF]. In this study, we will consider also leaf fall proxies because many autumn phenology studies use these proxies to associate a timing to the autumn senescence process at whole without considering the possible large time lag between leaf fall and leaf senescence onset [START_REF] Gill | Changes in autumn senescence in northern hemisphere deciduous trees: a metaanalysis of autumn phenology studies[END_REF]. Estimates of leaf senescence based on the seasonal trend of fallen leaf biomass or a defined percentage of fallen leaves (e.g. 50%) are common practice in the field because of their simple character [START_REF] Gill | Changes in autumn senescence in northern hemisphere deciduous trees: a metaanalysis of autumn phenology studies[END_REF]. For stands, proxies of the this fourth approach can be derived from remote sensing indices that can be related to leaf biomass, such as the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI; [START_REF] Rouse | Monitoring vegetation systems in the Great Plains with ERTS[END_REF][START_REF] Huete | A comparison of vegetation indices over a global set of TM images for EOS-MODIS[END_REF][START_REF] Wu | Contrasting responses of autumn-leaf senescence to daytime and night-time warming[END_REF].

Overall, it is important to note that remote sensing proxies for determination of leaf senescence onset in forests overcome sampling problems and issues related to visual observations. However, they also have inherent problems, not only related to their indirect character, spatial-and temporal resolution, but also to potential biases introduced by the variability of the vertical canopy profile during autumn [START_REF] Koike | Autumn coloring, photosynthetic performance and leaf development of deciduous broad-leaved trees in relation to forest succession[END_REF].

The accuracy, quality and comparability of the proxies described above to determine the leaf senescence onset and leaf fall timing are unknown. Moreover, the use of different proxies complicates the interpretation and comparison of previous results and meta-analyses [START_REF] Gill | Changes in autumn senescence in northern hemisphere deciduous trees: a metaanalysis of autumn phenology studies[END_REF]. There is an urgent need for a methodological comparison that could bring clarity in this field and put forward the most reliable and convenient approaches to estimate the leaf senescence onset in temperate deciduous forest trees.

In this study, we critically compared common proxies (both field-and remote sensing-based) of the approaches to measure the leaf senescence onset (leaf N content, chlorophyll degradation and leaf coloration) and leaf fall timing, for a European beech (Fagus sylvatica L.) stand at Kleine Schietveld (KS) close to Antwerp (Belgium) during late summer-autumn 2017 and 2018. Furthermore, field-and remote sensing-based proxies of chlorophyll degradation, leaf coloration and leaf fall were also applied in 2017 at a stand of silver birch (Betula pendula Roth.) at KS and at a stand of pedunculate oak (Quercus robur L.) at the nearby forest area Park of Brasschaat (PB). Finally, to evaluate the impact of different sites

Accepted Article

.

(KS and PB) on the detection of the leaf senescence onset, we compared for 2017 two field proxies (based on chlorophyll degradation and leaf coloration) for the abovementioned stands, plus one oak stand at KS, one beech stand at PB and one birch stand also at PB. As for several proxies the problems are not in the measurements but rather in the data analysis, for all proxies we compared a standard calculation procedure (based on threshold of 50% change in the seasonal variable used to derive the leaf senescence onset or leaf fall) with an improved method described in this study and based on breakpoint analyses (i.e. determination of shifting points in the seasonal pattern of a variable through piecewise linear regressions).

The general objectives of the study can be summarized by three questions: (1) what are the differences in estimates of leaf senescence onset and leaf fall timing resulting from the application of different proxies? (2) are these differences related to methodological issues and can they be amended? (3) what are the best applications for the different methods?

Materials and Methods

Study sites and species

The study was carried out at two locations close to Antwerp, Belgium: the Klein Schietveld in Kalmthout and Kapellen (51°21'N,4°37'E, 21 m.a.s.l) and the Park of Brasschaat in Brasschaat (51°12'N,"4°26'E, 10 m.a.s.l.). The PB is around eight km more southwards than the KS and its soil is moderately wetter and more nutrient rich than the soil at KS. Both locations experience a temperate maritime climate with an average annual temperature of 10.5 °C and an average rainfall of 919 mm yr -1 regularly distributed throughout the year [START_REF] Campioli | Current status and predicted impact of climate change on forest production and biogeochemistry in the temperate Accepted Article . oceanic European zone: review and prospects for Belgium as a case study[END_REF]. Five tropical days with a maximum temperature above 30 °C and six ice days with a maximum temperature below 0 °C are recorded on average annually. For both locations the wind mainly blows from South-West. The average amount of sun radiation and total hours of sunshine per year are 2.75 kWh m -2 d -1 and 1616 hours, respectively (KMI, 2010b;KMI, 2010a).

In 2017, the weather was relatively normal with a warm summer (June to August) and a normal autumn (September to November). In summer, the average temperature was 18.6 °C and the total rainfall was 179.9 mm (KMI, 2017b). In autumn, the average temperature was 11.3 °C and the total rainfall was 226.5 mm (KMI, 2017a). Unlike 2017, the weather in 2018 was exceptional with the occurrence of heat waves and drought events during summer that either broke local weather records or that only reoccur every 30 years. The summer was

Accepted Article

. exceptionally warm and dry as the average temperature was 19.8 °C and the total rainfall of 134.7 mm fell during only 20 rainy-days (KMI, 2018b). While the average autumn temperature was similar to those in 2017 (11.8°C), the 2018 autumn was also abnormally dry and the total rainfall of 168.5 mm fell only during 32 rainy-days (KMI, 2018a).

For this study, three tree species were selected: beech, oak and birch. Oak and particularly beech have a closed canopy and a deterministic growth pattern (i.e. one or two springsummer leaf flushes) with senescence starting first in their outer sun-exposed leaves [START_REF] Koike | Autumn coloring, photosynthetic performance and leaf development of deciduous broad-leaved trees in relation to forest succession[END_REF]). Birch has an open canopy and a non-deterministic growth pattern (i.e. continuous leaf flushing) with leaf senescence starting in the older leaves in the inner parts of the canopy [START_REF] Koike | Autumn coloring, photosynthetic performance and leaf development of deciduous broad-leaved trees in relation to forest succession[END_REF]. Because of the different degree of canopy openness, we classified beech and oak leaves as sun-and shade-leaves, whereas birch leaves were always considered sunleaves.

Our main stand (beech KS) has an extension of 3-ha (ca. 150 × 200 m) and is monospecific and homogeneous. The two other intensively study stands (birch KS and oak PB, 0.8 ha and 3 ha, respectively) are also monospecific and homogeneous. On the other hand, the three ancillary stands, are smaller (<0.5 ha as beech PB) or with mixed species (as birch PB and oak KS). The beech and oak trees in the study sites were mostly planted but little human management has been performed, particularly in the last decades. Beech trees were 60 to 70 years old and oaks were 60 to 120 years old. As for the birch sites, the trees were naturally established and were 50 to 60 years old. At each stand, eight individuals were selected for canopy measurements and four for leaf traits measurements, except for the beech trees at KS, where the amount of trees selected was double (thus sixteen and eight, respectively). All individual trees were selected for their dominance and vitality.

Field-based proxies

Referring to the four general approaches (leaf nutrient relocation, chlorophyll degradation, leaf coloration and leaf fall), the investigated field-based proxies are: the leaf N content (first approach); the chlorophyll content index (CCI, second approach); loss of canopy greenness (third approach) and fallen leaf biomass (fourth approach). An overview of the proxies is given in Table 2.

Accepted Article

.

Leaf traits: Chlorophyll and N content

Leaf traits were measured from the end of July to the end of November 2017 and 2018.

Leaves from each tree were collected by tree-climbers on eight occasions: once every two to three weeks throughout the study period. Measurement days were generally dry and sunny, and sampling took place from 10 a.m. to 4 p.m. For beech and oak, five sun-leaves and five shade-leaves were collected. For birch, five sun-leaves were collected. The chlorophyll content index of these leaves was immediately measured after collection with a chlorophyll content meter (CCM-200 plus, OPTI-SCIENCES Inc., USA). All collected leaves were brought to the University of Antwerp. After removing the leaf petiole, the leaves were airdried at 70 °C for 48 hours and individually ground to pass a 0.5-mm sieve in an ultracentrifugal mill (Model ZM 200, Retsch GmbH, Haan, Germany). The resulting fine powder was used to determine the total N content by dry combustion based on the Dumas method using an elemental analyzer (Model FLASH 2000, Thermo Fisher Scientific, Waltham, MA, USA). Because of the small available biomass, sun-and shade-leaves were pooled in the latter analysis (which was done only for beech trees at KS).

To compare the CCI trend with actual spectrophotometrically measured chlorophyll values, two circular samples of leaf tissue ("circles") were collected on one occasion per month at beech KS from each leaf sampled for CCI using a shaped cylinder of 10 mm diameter. For each leaf, (i) one circle was weighed after 48 hours at 70°C to determine specific leaf area, while (ii) the second circle (stored before processing at -80°C) was ground in a centrifuge with glass beads and dissolved in ethanol. Subsequently, the absorption of the resulting supernatant was measured with a spectrophotometer (Smart Spec Plus Spectrophotometer, Bio-Rad Laboratories, USA) at wavelengths of 662 nm and 644 nm for chlorophyll a and chlorophyll b, respectively. Chlorophyll concentrations where obtained using the following formula [START_REF] Holm | Chlorophyll Mutations in Barley[END_REF][START_REF] Vonwettstein | Chlorophyll-letale und der submikroskopische Formwechsel der Plastiden[END_REF]: eq. 1 with E662 and E644 being the measured absorption value of chlorophyll a and b, respectively. The CCI measurements showed to have a linear relationship with the chlorophyll concentration measurements (R²=0.58, Fig. S1) and to capture its seasonal pattern (Fig. S2).

Canopy characteristics: Loss of canopy greenness and fallen leaf biomass

Every week from early September till late November 2017 and 2018 (date; t), the loss of canopy greenness (xt) of each of the eight to sixteen trees was visually estimated following

Accepted Article

. the method of [START_REF] Vitasse | Responses of canopy duration to temperature changes in four temperate tree species: relative contributions of spring and autumn leaf phenology[END_REF]. In late summer and early autumn xt was estimated directly as a percentage. Later in the season, when the process became more intense, xt was estimated through a combined rating of the percentage of leaves that had changed color (αt) and the percentage of leaves that had fallen (βt; fallen leaves were assumed to be colored).

They were related to xt with the following formula:

eq. 2

The fallen leaf biomass was estimated for eight beech trees at KS. Eight rectangular litter traps of ca. 0.2 m 2 (53 × 37 cm) were placed aboveground around each individual tree (Fig. S3). They were placed pairwise in each direction of the wind: for each pair of traps, one was put three meters away from the tree stem, whereas the other was put 1.5 m away from the stem. The litter traps were emptied every week from early October to late November, and their litter biomass was measured after air-drying at 70 °C in an oven for 48 hours. The data of the leaf traps were used for the leaf litter proxy.

Remote sensing indices

Remote sensing proxies were derived from European Space Agency Sentinel-2 TOC (Top Of Canopy) images with ten or twenty meter pixel resolution depending on the spectral bands.

First, the coordinates of each tree were measured (with R8s GNSS receiver and TSC3 controller, Trimble Inc., USA). Second, for each tree, the coordinates were associated with the corresponding Sentinel 2 pixel. Third, remote sensing indices were derived for this representing the five used Sentinel-2 bands, were calculated using the following formulas [START_REF] Hawryło | Estimating defoliation of Scots pine stands using machine learning methods and vegetation indices of Sentinel-2[END_REF]: eq. 3

Accepted Article

. eq. 4 eq. 5 eq. 6 eq. 7

The Sentinel-2 bands cover the following part of the wavelength spectrum; B02: 459.1 to 525.1 nm, B04: 649.4 to 680.5 nm; B05: 695.8 to 711.8 nm; B06: 731.6 to 746.6 nm and B08: 779.9 to 885.9 nm. B02, B04 and B08 have a ten meter spatial resolution, while B05 and B06 have a twenty meter resolution.

The CHL-RED-EDGE and the MTCI use the spectral bands related to the vegetation canopy chlorophyll absorption feature. They were both proposed and validated as robust, sensitive and easy-to-measure indices for chlorophyll content in vegetation [START_REF] Gitelson | Remote estimation of canopy chlorophyll content in crops[END_REF][START_REF] Dash | Evaluation of the MERIS terrestrial chlorophyll index (MTCI)[END_REF][START_REF] Dash | The MERIS terrestrial chlorophyll index[END_REF][START_REF] Clevers | Remote estimation of crop and grass chlorophyll and nitrogen content using red-edge bands on Sentinel-2 and -3[END_REF]. The PSRI has been suggested as an accurate leaf senescence proxy, sensitive to carotenoids in leaves (Merzlyak et al., 1999;[START_REF] Rautiainen | Seasonal reflectance dynamics of common understory types in a northern European boreal forest[END_REF][START_REF] Cole | Spectral monitoring of moorland plant phenology to identify a temporal window for hyperspectral remote sensing of peatland[END_REF][START_REF] Ren | Assessing plant senescence reflectance index-retrieved vegetation phenology and its spatiotemporal response to climate change in the Inner Mongolian Grassland[END_REF]. While the NDVI and the EVI are similar indices using near-infrared and visible wavelengths for estimation of the green biomass, the EVI uses additional reflected wavelengths to correct for the effect of the atmosphere, soil conditions or the solar incidence angle that can influence the reflectance values [START_REF] Rouse | Monitoring vegetation systems in the Great Plains with ERTS[END_REF][START_REF] Zhou | Variations in northern vegetation activity inferred from satellite data of vegetation index during 1981 to 1999[END_REF][START_REF] Jiang | Development of a two-band enhanced vegetation index without a blue band[END_REF][START_REF] Zhang | Monitoring fall foliage coloration dynamics using time-series satellite data[END_REF].

Data analysis: determination of leaf senescence onset and leaf fall

Following the classical approach, for all proxies, the date of leaf senescence onset and leaf fall timing was determined as the day in which the seasonal value of the different variables reached a fixed threshold of 50% of the summer maximum [START_REF] White | Intercomparison, interpretation, and assessment of spring phenology in North America estimated from remote sensing for 1982-2006[END_REF]. On the other hand, to improve these estimates, we calculated also the leaf senescence onset and leaf fall timing as the day when the variable used for the assessment started to decline, or increase, substantially in early autumn. This was done through the use of a 'breakpoint analysis' (i.e. piecewise regression; [START_REF] Rosenthal | Photosynthetic decline and pigment loss during autumn foliar senescence in western larch (Larix occidentalis)[END_REF][START_REF] Galvagno | Seasonal course of photosynthetic efficiency in Larix decidua Mill. in response to temperature and change in pigment composition during senescence[END_REF][START_REF] Menzel | Patterns of late spring frost leaf damage and recovery in a European beech (Fagus sylvatica L.) stand in south-eastern Germany based on repeated digital photographs[END_REF][START_REF] Wingate | Interpreting canopy development and physiology using a European phenology camera network at flux sites[END_REF]. The breakpoint correspond to the date when the seasonal course of the

Accepted Article

. measured variable can be 'broken' into two linear regressions with different slopes or even sign (Fig. S4). This method is coded in the R-package 'segmented', a package to fit regressions with broken-line relationships [START_REF] Vito | segmented: an R Package to Fit Regression Models with Broken-Line Relationships[END_REF][START_REF] Wickham | ggplot2: Elegant Graphics for Data Analysis[END_REF][START_REF] Wickham | dplyr: A Grammar of Data Manipulation[END_REF].

Statistical analyses

All statistical analyses were performed with the program R, version 3.5.0 (R Core Team, 2019). Graphs were built using the packages 'dplyr' and 'ggplot2'. All response variables were tested for normality following [START_REF] Zuur | A protocol for conducting and presenting results of regression-type analyses[END_REF]. For all analyses, the uncertainty reported here refers to inter-tree variability.

Analysis to compare the different proxies. This analysis was done on data of 2017 and 2018, separately, of beech KS, considering both field-and remote sensing-based proxies and only using breakpoint analysis as calculation method. The requirement of homoscedasticity was not met. Thus, to test for differences in the date of leaf senescence onset among all different proxies, Welch's analysis of variance test was performed with the measurement proxy as predictor variable and the date of leaf senescence onset (or leaf fall) as response variable. To test the level of significance among the proxies, Dunn's post-hoc test was performed using the R packages 'dunn.test' and 'FSA' [START_REF] Dinno | test: Dunn's Test of Multiple Comparisons Using Rank Sums[END_REF][START_REF] Ogle | R: A language and environment for statistical computing[END_REF]. Afterwards, a Bonferonni correction was applied to adjust for multiple comparisons.

Analysis to compare the different calculation methods (breakpoint vs 50% threshold). This analysis was performed on combined data from 2017 and 2018 for beech KS. To test for the impact of the calculation method on the leaf senescence onset and leaf fall timing, we constructed a linear model using as predictor variable the calculation method and as response variable the date of leaf senescence onset or leaf fall timing. Since the requirement of homoscedasticity was not met, Welch's analysis of variance test was used. For each proxy, to test the level of significant difference between the calculation methods, t-tests were used.

Analysis to compare the different species and sites. First, we wanted to check whether the potential differences in the date of leaf senescence onset and leaf fall obtained by using the different proxies depended on the species. Therefore, the analysis on the different proxies done for beech KS (see above) was repeated for birch KS and oak PB using data of 2017.

Second, we wanted to verify if site (PB vs KS) had an impact on the results provided by the different proxies. This analysis was performed for all stands in KS and PB, considering only the proxies available for all six stands, i.e. the field-based proxies of CCI and loss of canopy Accepted Article . greenness. In the latter analysis, the requirements of normality and homoscedasticity were not met. Thus, to test for site or species effect on the date of leaf senescence onset per proxy, we made a linear model with as predictor variables the proxy (CCI, loss of canopy greenness), species, site, the interaction between the proxy and the species and the interaction between the proxy and the site. The response variable was the date of leaf senescence onset. We note that 'site differences' comprise not only differences in site fertility and environmental conditions but also differences in tree age, stand structure etc.

Results

Comparison among proxies and years.

The 2017 analysis shows no significant differences among the results provided by the proxies of leaf senescence onset (Fig. 1 panel B; Fig. 2 panels A, B, D, E andG; Table S1). The mean values of leaf senescence onset that these proxies provided were between the 6 th of October (DOY = 279 ± 7) and the 22 th of October (DOY = 295±5). The mean onset date of the litter proxy (DOY =315±1; 11 th of November) was only significantly later than the EVI and PSRI (Fig. 1; Fig. 2 panels C, F and I; Table S1). The NDVI provided an intermediate value of leaf senescence onset (DOY = 300±7) which was not significantly different compared to the results from the other proxies of leaf senescence onset and leaf fall timing. However, if pairwise comparisons (thus without the need of considering Bonferroni correction) were done between results from litter proxy and from one of the leaf senescence onset proxies, the differences were consistently significant (Table S1).

The results of the analysis performed on data of 2018 showed a similar pattern as the results from the data of 2017 (Fig. 1 panel B andD). However, a major exception was represented by the results obtained with the loss of canopy greenness proxy. In fact, the latter proxy provided a significantly earlier date of leaf senescence onset (DOY = 248±7) than the other proxies, in particular, the CCI and CHL-RED-EDGE (Fig. 1 panel D; Fig. 2 panel G). Onset of leaf senescence measured through the loss of canopy greenness was 40 days earlier in 2018 than 2017 (t.test; p < 0.05). On the other hand, the value of leaf senescence onset obtained with the CCI proxy in 2018 was surprisingly similar to those of 2017 (t.test; p = 0.52). As for 2017, significant differences were expressed more when Bonferroni correction was not considered (see Table S1).

Accepted Article

.

Comparison between calculation methods

Considering the data of 2017 and 2018 for the beech trees at KS, the leaf senescence onset and leaf fall timing were differently estimated by the 50% threshold method and the breakpoint analysis (Welch's ANOVA; p < 0.001). The extent of this difference depended on the proxy (Fig. 3; Table S2). In seven cases out of nine (i.e. all proxies except for the NDVI and litter proxies), the mean leaf senescence onset derived from the 50% threshold analysis and the breakpoint analysis were significantly different (t.test; p < 0.05), with the 50% threshold analysis providing later estimates, ranging from nine days to five weeks.

Comparison among species and sites

Birch presented a very similar pattern as for beech, with no differences in results from the different proxies in 2017 (Fig. 1 panels A andB; Table S1 andS3, and Fig. S5). On the other hand, at the oak stand, results from the EVI were significantly different than results from the field-based methods and the NDVI provided later dates than the EVI, MTCI and PSRI (Fig. 1 panel C, Fig. S6). However, across species: (i) the field-based proxies CCI and loss of canopy greenness did not provide different results between each other and (ii) the CHL-RED-EDGE and the NDVI provided results consistently similar to the ones obtained with the field-based proxies. Corroborating the former point, data from the CCI and loss of canopy greenness for all six stands showed that the interaction between proxy and species, and proxy and site was not significant. Thus, site and species did not affect the results of the proxies.

Discussion

We presented an extensive comparison among different proxies to measure the leaf senescence onset and leaf fall timing in temperate deciduous forest trees. We did not only compare the proxies based on different (though interconnected) processes (i.e. nutrient relocation, chlorophyll degradation, leaf coloration and leaf loss) but we also compared fieldversus remote sensing proxies. We also presented the use of the breakpoint analysis to determine the leaf senescence onset. This breakpoint analysis is useful when the seasonal change (decline or increase) of the variable under consideration (e.g. chlorophyll content, loss of canopy greenness) is not consistently clear due its variability. The use of the breakpoint analysis was extended to remote sensing indices and showed successful results.

Considering the field-based proxies, the two most direct methodologies to detect the leaf senescence onset (i.e. the approaches based on leaf N relocation and on chlorophyll degradation), provided similar results even when using indirect measurements of chlorophyll (e.g. CCI). Therefore, the choice between these two approaches depends on the researcher

Accepted Article

. preference and practical aspects. Our data of 2017 showed that the proxy based on leaf coloration could be used as valid alternative to the chlorophyll degradation proxy because results from both approaches were not significantly different. However, our data of 2018 showed otherwise, which is surprising since the onset of leaf senescence from the CCI proxy were similar in both years. The likely explanation is that in 2018 many leaves were damaged or fallen earlier than usual due to the severe drought and extreme temperatures in early summer. This process is called 'accelerated leaf senescence' and it is substantially different than the autumn senescence studied here as it does not involve leaf nutrient relocation, but just leaf damage and death [START_REF] Gunthardt-Goerg | Linking stress with macroscopic and microscopic leaf response in trees: new diagnostic perspectives[END_REF][START_REF] Estiarte | Alteration of the phenology of leaf senescence and fall in winter deciduous species by climate change: effects on nutrient proficiency[END_REF]. The leaves that were measured with the chlorophyll content meter in 2018 were likely healthy leaves which followed the normal cues for onset of the autumn senescence process.

Therefore, it appears that the loss of canopy greenness proxy should not be used when trees are subjected to severe stress (e.g. extreme drought, heat haves, pests) which cause substantial accelerated leaf senescence. On the other hand, both proxies (CCI and loss of canopy greenness) were not species-or site-dependent, which is relevant for their application.

It is important to stress that the agreement among the three general field-based approaches to detect leaf senescence onset (leaf N relocation, chlorophyll degradation and leaf coloration) was obtained when applying the new calculation method based on breakpoints. In fact, the dates of leaf senescence onset calculated with the 50% threshold -which is the more established calculation method-were around one-and-a-half to five weeks later than those based on the breakpoint analysis and with significant differences between proxies. Finally, our study showed that field proxies based on leaf litter introduce a systematic time lag compared to the proxies of leaf senescence onset (three to four weeks) and remote sensing proxies related to leaf fall such as EVI and NDVI (two to four weeks). Quantifying this time lag is relevant as, in many studies, leaf fall timing is taken as an indication of the whole leaf senescence process [START_REF] Gallinat | Autumn, the neglected season in climate change research[END_REF][START_REF] Gill | Changes in autumn senescence in northern hemisphere deciduous trees: a metaanalysis of autumn phenology studies[END_REF]. Also, it should be taken into account that leaf fall can occur in early in the season (e.g. due to damage by high temperatures, droughts, winds, see above).

Our study showed that remote sensing proxies can be used to detect the leaf senescence onset with accuracy, in particular the CHL-RED-EDGE and NDVI. The CHL-RED-EDGE is known to be a reliable index related to chlorophyll degradation [START_REF] Gitelson | Remote estimation of canopy chlorophyll content in crops[END_REF][START_REF] Clevers | Remote estimation of crop and grass chlorophyll and nitrogen content using red-edge bands on Sentinel-2 and -3[END_REF]. However, until now, it was uncertain whether it was a reliable proxy for the detection of the leaf senescence process. Despite the not univocal connection between the

Accepted Article

. chlorophyll degradation and the process of nutrient relocation, which characterizes the autumn leaf senescence process, our study shows that the CHL-RED-EDGE can indeed be a reliable proxy for the leaf senescence onset. The EVI and NDVI have both been tested for monitoring autumn phenology extensively and found to be closely related to seasonal changes in leaf biomass, which should provide estimates of leaf senescence onset with systematic delay (see above; [START_REF] Wang | On the relationship of NDVI with leaf area index in a deciduous forest site[END_REF][START_REF] Croft | Temporal disparity in leaf chlorophyll content and leaf area index across a growing season in a temperate deciduous forest[END_REF]Klosterman, Stephen & Richardson, Andrew D., 2017). However, this was not the case as, especially the NDVI, provided consistent results with the field proxies of leaf senescence onset. Note that annual curves of both NDVI and EVI reflect also leaf coloration, which would explain the match with field proxies of leaf senescence onset [START_REF] Yang | Beyond leaf color: Comparing camera-based phenological metrics with leaf biochemical, biophysical, and spectral properties throughout the growing season of a temperate deciduous forest[END_REF]. On the other hand, the EVI did not compare well with field methods for oak. Maybe the corrections, which are present in the EVI for errors related to vegetation density and its better correlation to chlorophyll changes, made it too sensitive to coloration changes in this species [START_REF] Zhou | Variations in northern vegetation activity inferred from satellite data of vegetation index during 1981 to 1999[END_REF][START_REF] Zhang | Monitoring fall foliage coloration dynamics using time-series satellite data[END_REF].

It is not possible to define the time resolution needed to obtain the best results with the breakpoint analysis, as this analysis is also affected by the data quality of the time-series, their seasonal trend, the presence of outliers etc. However, because the breakpoint analysis is particularly problematic in case of abrupt changes, we recommend a weekly (or fortnightly) time resolution with a start and end of the observations at least three to four weeks (or oneand-a-half to two months in case of fortnightly assessment) before and after the expected breakpoint date.

In this study we did not test the performance of the proxy based on the coloration change measured with repeated colour images of the forest canopy throughout the season [START_REF] Cai | Quantifying the Onset and Progression of Plant Senescence by Color Image Analysis for High Throughput Applications[END_REF]Klosterman, S. & Richardson, A. D., 2017). This method can offer high temporal resolution and reduces the subjective character of the classical visual coloration observations, but it requires the installation of a permanent camera on a tall mast or frequent drone measurements, both options considered unpractical at our sites [START_REF] Chipkin | Effect of propranolol on antinociceptive and withdrawal characteristics of morphine[END_REF][START_REF] Brown | Using phenocams to monitor our changing Earth: toward a global phenocam network[END_REF][START_REF] Liu | Using canopy greenness index to identify leaf ecophysiological traits during the foliar senescence in an oak forest[END_REF].

Conclusion

We provide here an answer to our three general research questions. ( 1) What are the differences in onset of leaf senescence resulting from the application of different proxies?

Classical methods based on 50% threshold in the seasonal value of the variable used to derive

Accepted Article

. leaf senescence onset introduce a significant lag in the detection of leaf senescence.

Therefore, these methods should not be used. Instead, the application of breakpoint analysis provides more reliable leaf senescence estimates and no difference among approaches based on leaf N relocation, chlorophyll degradation and leaf coloration. However, the proxy based on loss of canopy greenness should not be used under severe stress conditions damaging significantly the canopy. The leaf fall date from data of the leaf litter is three to four weeks later than leaf senescence onset. (2) Are differences among methods related to methodological issues and can they be amended? The use of thresholds for the seasonal value of the variable used to derive leaf senescence should be substituted by other calculation methods that better consider the seasonal changes leading to leaf senescence onset or leaf fall.

We proposed here the breakpoint analysis as alternative, but other approaches (e.g. derivatives to determine peaks and trend changes, spline interpolation or Bayesian and Pruned Exact Linear Time change point techniques) are also possible depending on the time resolution and trend (i.e. gradual or abrupt) of the available data [START_REF] Zhang | Monitoring vegetation phenology using MODIS[END_REF][START_REF] Klosterman | Evaluating remote sensing of deciduous forest phenology at multiple spatial scales using PhenoCam imagery[END_REF][START_REF] Liu | Using canopy greenness index to identify leaf ecophysiological traits during the foliar senescence in an oak forest[END_REF][START_REF] Richardson | Tracking vegetation phenology across diverse North American biomes using PhenoCam imagery[END_REF]. ( 3) What are the best applications for the different methods? (i) All three key field-based proxies, leaf N content, CCI and loss of canopy greenness could be used for the most fine-scale and detailed applications of leaf senescence onset data (e.g. physiological studies investigating single tree phenology with molecular methods and eco-physiological studies relating phenology to tree growth) in absence of severe stress. However, we note that the rarely used proxy based on leaf N content was the one providing results with the lowest intra-individual variability, which might be advantageous in these types of studies. (ii) CCI is more suited than the loss of canopy greenness for ecological application and long term monitoring of single trees, with the caveat of the practical problems to access the leaves. (iii) The CHL-RED-EDGE and NDVI could be reliably used to determine leaf senescence onset at stand level (thus for large scale applications) across species and stress levels. (iv) Data from field-based proxies of leaf fall dynamics might be used in factorial studies, when also the control is measured in this way. Meta-analyses on historical time-series of leaf fall data should take into account that these data refer to the end of the leaf senescence process and are sensitive to random external factors (e.g. strong winds). Overall, our results improve the setting up of autumn phenology research, the elucidation of previous contrasting findings and the development of more accurate monitoring protocols of forest ecosystems. Table S1: Differences between results from pair of proxies of leaf senescence onset or leaf fall timing applied to the beech stand of Klein Schietveld.

Brief legends for Supporting Information

Table S2: Significance level of the difference between data of leaf senescence onset (or leaf fall timing) determined with two calculation methods i.e. breakpoint analysis and 50% threshold method, for nine assessment proxies applied at the beech stand of Klein Schietveld.

Table S3: Differences between results from pair of proxies of leaf senescence onset or leaf fall timing applied separately to the birch stand of Klein Schietveld and the oak stand at Park of Brasschaat.

Accepted Article

. for Fagus sylvatica at Klein Schietveld (Antwerp, Belgium) determined with two calculation methods i.e. breakpoint (blue) and the 50% threshold method (red; see the 'Data analysis: determination of leaf senescence onset and leaf fall' subsection), and eight assessment proxies. The dates were averaged over two years (2017 and 2018), except for the leaf N proxy for which only data of 2017 was available. Dots indicate that the respective proxy refers to the leaf senescence onset, while diamond shapes indicate that the respective proxy refers to the timing of the leaf fall. The notation between pairs of values for the same proxy indicates a significant difference (***: p < 0.001; **: p < 0.01; *: p < 0.05) or no significant difference (n.s.) checked with t-test. The black horizontal lines were added to easily detect the pair of methods for each proxy. MTCI, MERIS terrestrial reflectance index; EVI, enhanced vegetation index; CHL-RED-EDGE, chlorophyll red-edge index; CCI, chlorophyll content index; PSRI, plant senescence reflectance index; NDVI, normalized difference vegetation index.

Accepted Article

.

Accepted Article

.

  (central) pixel and the eight surrounding pixels. Finally, the values of the nine pixels were averaged and used as such in this analysis. The indices were calculated for the period of the field sampling, from the end of July to late November 2017 and 2018. On average, Sentinel-2 observations were available at a five day frequency. Five indices were extracted: the chlorophyll red edge index (CHL-RED-EDGE) and the MERIS terrestrial chlorophyll index (MTCI) that indirectly correlate with chlorophyll degradation; the plant senescence reflectance index (PSRI), indirectly correlated with canopy coloration, and the normalized difference vegetation index (NDVI), and the enhanced vegetation index (EVI), indirectly correlated with leaf fall. The remote sensing indices, with B02, B04, B05, B06 and B08

Figure S1 :

 S1 Figure S1: Relationship between chlorophyll content index and chlorophyll concentration of beech leaves.

Figure S2 :

 S2 Figure S2: Comparison between the seasonal trend of the chlorophyll content index and chlorophyll concentration from spectroscopical analysis of beech leaves.

Figure S3 :

 S3 Figure S3: Schematic representation of the litter trap placement at the Fagus sylvatica stand at Klein Schietveld.

Figure S4 :

 S4 Figure S4: Schematic example of how the leaf senescence onset was calculated through the breakpoint and the 50% threshold method performed on leaf nitrogen content data of all Fagus sylvatica trees at Klein Schietveld.

Figure S5 :

 S5 Figure S5: Seasonal pattern of the variables used to derive leaf senescence onset and leaf fall timing for Betula pendula.

Figure S6 :

 S6 Figure S6: Seasonal pattern of the variables used to derive leaf senescence onset and leaf fall timing for Quercus robur.

FIGURESFigure 1 :

 1 FIGURESFigure1: Dates of leaf senescence onset and leaf fall timing (means with error bars representing ± 1 standard error; n = 4-16; the standard error shows the inter-individual variability) for Fagus sylvatica and Betula pendula trees at Klein Schietveld, and Quercus robur trees at Park of Brasschaat (both close to Antwerp, Belgium) calculated with breakpoint analyses. For Betula (panel A) and Quercus (panel C), leaf senescence onset and leaf fall timing were calculated for five and two assessment methods, respectively, for 2017. For Fagus, leaf senescence onset and leaf fall timing were calculated for six and three assessment methods, respectively, in 2017 (panel B) and 2018 (panel D). Red and blue shapes indicate field-and remote sensing proxies, respectively. Dots indicate that the respective proxy refers to the leaf senescence onset, while diamond shapes indicate that the respective proxy refers to the timing of the leaf fall. Different capital letters on the right side indicate significant differences among proxies by a Dunn's test (p < 0.05) after Bonferroni correction. DOY, day of the year; MTCI, MERIS terrestrial reflectance index; EVI, enhanced vegetation index; CHL-RED-EDGE, chlorophyll red-edge index; CCI, chlorophyll content index; PSRI, plant senescence reflectance index; NDVI, normalized difference vegetation index.

Figure 2 :

 2 Figure 2: Seasonal pattern of the variables used to derive leaf senescence onset and leaf fall timing for Fagus sylvatica in 2017 (red) and 2018 (blue). For each panel, dots represent means, error bars represent the standard error showing the inter-individual variability, the solid vertical line represents the mean date of leaf senescence onset or leaf fall timing, while the dotted horizontal lines represent its standard errors. MTCI, MERIS terrestrial chlorophyll index; EVI, enhances vegetation index; CCI, chlorophyll content index; CHL-RED-EDGE, chlorophyll red-edge index; PSRI, plant senescence reflectance index; NDVI, normalized difference vegetation index.

Figure 3 :

 3 Figure3: Dates of leaf senescence onset and leaf fall timing (means with error bars representing ± 1 standard error; n = 4-16; the standard error shows the inter-individual variability) for Fagus sylvatica at Klein Schietveld (Antwerp, Belgium) determined with two calculation methods i.e. breakpoint (blue) and the 50% threshold method (red; see the 'Data analysis: determination of leaf senescence onset and leaf fall' subsection), and eight assessment proxies. The dates were averaged over two years(2017 and 2018), except for the leaf N proxy for which only data of 2017 was available. Dots indicate that the respective proxy refers to the leaf senescence onset, while diamond shapes indicate that the respective proxy refers to the timing of the leaf fall. The notation between pairs of values for the same proxy indicates a significant difference (***: p < 0.001; **: p < 0.01; *: p < 0.05) or no significant difference (n.s.) checked with t-test. The black horizontal lines were added to easily detect the pair of methods for each proxy. MTCI, MERIS terrestrial reflectance index; EVI, enhanced vegetation index; CHL-RED-EDGE, chlorophyll red-edge index; CCI, chlorophyll content index; PSRI, plant senescence reflectance index; NDVI, normalized difference vegetation index.

  FIGURESFigure1

  Figure 2

Figure 3

 3 Figure 3

Table 1 :

 1 Overview of the main three approaches to determine leaf senescence onset in deciduous forests.

	TABLES
	Approach	Principle/process of reference	Advantages	Disadvantages
	1) Leaf N content	Nutrient relocation from	-Direct method	-Decline of leaf N in autumn
		leaves to perennial organs		not consistently clear
				-Laborious
				-No direct remote sensing
				proxies
	2) Chlorophyll degradation Decline of chlorophyll in	-Direct method	-Decline of chlorophyll in
		Autumn	-Remote sensing approaches possible	autumn not consistently clear
				-Laborious
	3) Leaf coloration	Measurement of coloration	-Simple	-Possibly subjective
		appearing as consequence	-Remote sensing approaches possible	-Possible systematic lag
		of chlorophyll degradation

Table 2 :

 2 Characteristics of the nine proxies of leaf senescence onset and leaf fall timing compared in this study.On average between late July and Late November 2017 and 2018 at the study sites.

	Proxy	Abbreviation	Process	Time resolution	Spatial resolution
	Leaf senescence onset
	Nitrogen	N	Nutrient relocation	Monthly		tree
	Chlorophyll content index	CCI	Chlorophyll decline	fortnightly		tree
	Loss of canopy greenness	-	Leaf coloration and fall	Weekly		Tree
	MERIS terrestrial chlorophyll index	MTCI	Chlorophyll decline	Every five days	a	Stand
	Chlorophyll red edge index	CHL-RED-EDGE	Chlorophyll decline	Every five days	a	Stand
	Plant senescence reflectance index	PSRI	Leaf coloration	Every five days	a	Stand
	Leaf fall
	Normalized difference vegetation index	NDVI	Biomass decline	Every five days	a	Stand
	Enhanced vegetation index	EVI	Biomass decline	Every five days	a	Stand
	Litter	-	Leaf fall	Weekly		Tree -stand
	a

Acknowledgments

This research was funded by the ERC Starting Grant LEAF-FALL (714916) and the University of Antwerp (DOCPRO4 'Determination of the drivers of the onset of autumn leaf senescence in temperate deciduous forests: the relationship between leaf dynamics, tree growth and photoperiod'). Campioli M. was a Postdoctoral Fellow of the Research Foundation-Flanders (FWO) and Balzarolo M. acknowledges the support provided by the EU Horizon 2020 Research and Innovation program under the Marie Skłodowska-Curie grant (INDRO, grant no. 702717). AbdElgawad H. was supported by a postdoctoral fellowship from the FWO (grant no. 12U8918N). Further funding was given by the Methusalem funding scheme of the Flemish Community through the Research Council of the University of Antwerp and by the Flemish Science Foundation (FWO, Brussels). We would like to thank the Belgian institutions that gave permission to conduct research in the study areas: Agency for Forest and Nature of the Flemish Government (ANB), the Military Defense of Belgium (Defensie), City of Brasschaat. Special thanks are due to Dirk Leyssens (ANB).

Author Contributions

All authors contributed to discussions and revisions.