Kaaviya Velumani
email: kvelumani@hiphen-plant.com

Simon Madec

Benoit De Solan

Raul Lopez-Lozano

Jocelyn Gillet

Jeremy Labrosse

Stephane Jezequel

Alexis Comar

Frédéric Baret

An automatic method based on daily in situ images and deep learning to date wheat heading stage

Keywords: Phenology, Internet of things for Agriculture, Convolutional Neural Networks, field sensors, phenology modelling

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

Introduction

Phenological observations are essential in agronomy, as crop management strategies (irrigation, fertilizing or crop protection) are planned considering plant development [START_REF] Brown | Nitrogen management for hard wheat protien enhancement[END_REF][START_REF] Chmielewski | Phenology in Agriculture and Horticulture[END_REF]. In wheat, the heading stage is one of the critical developmental phases as the plant becomes highly sensitive to abiotic stress -heat stress, frost, water constraints -with a strong impact on yield components [START_REF] Slafer | Sensitivity of Wheat Phasic Development to Major Environmental Factors: a Re-Examination of Some Assumptions Made by Physiologists and Modellers[END_REF]. Several studies have quantified the effect of post-heading abiotic stress on yield [START_REF] Ferris | Effect of High Temperature Stress at Anthesis on Grain Yield and Biomass of Field-grown Crops of Wheat[END_REF][START_REF] Gooding | Effects of Restricted Water Availability and Increased Temperature on the Grain Filling, Drying and Quality of Winter Wheat[END_REF][START_REF] Wheeler | The duration and rate of grain growth, and harvest index, of wheat ( Triticum aestivum L.) in response to temperature and CO 2[END_REF] due to a significant decrease in the grain weight,, along with a fall in the number of grains per plant. More recently, [START_REF] Balla | Heat stress responses in a large set of winter wheat cultivars (Triticum aestivum L.) depend on the timing and duration of stress[END_REF] analysed the possible impact of temperature on yield components at different development stages over a large wheat genotype panel, highlighting the high sensitivity of grain number to heat stress episodes around heading. Also, [START_REF] García | High night temperatures during grain number determination reduce wheat and barley grain yield: a field study[END_REF] and [START_REF] Lobell | Impacts of day versus night temperatures on spring wheat yields: A comparison of empirical and CERES model predictions in three locations[END_REF] have shown that the increase in minimum temperature during nights owing to changing climatic conditions accelerates the rate of crop development and reduces yield.

The timing and duration of stress through wheat developmental phases are thus essential to understand the impact of environmental factors on yield [START_REF] Sadras | Environmental modulation of yield components in cereals: Heritabilities reveal a hierarchy of phenotypic plasticities[END_REF]. Wheat phenology is driven by several eco-physiological mechanisms involving the response to temperature, photoperiod and vernalization [START_REF] Gate | Ecophysiologie du ble[END_REF] that are regulated by complex genetic pathways [START_REF] Guedira | Heading Date QTL in Winter Wheat (Triticum aestivum L.) Coincide with Major Developmental Genes VERNALIZATION1 and PHOTOPERIOD1[END_REF][START_REF] Whittal | Allelic variation of vernalization and photoperiod response genes in a diverse set of North American high latitude winter wheat genotypes[END_REF]. Phenology has been traditionally one of the most important traits used to genetically improve wheat adaptation, matching crop development -particularly reproductive and grain-filling phases -to the optimal growing conditions of a target environment [START_REF] Foulkes | Raising yield potential of wheat. III. Optimizing partitioning to grain while maintaining lodging resistance[END_REF][START_REF] Slafer | Wheat development: its role in phenotyping and improving crop adaptation[END_REF]. Indeed, phenology constitutes, in wheat breeding programs, one of the main levers enabling to optimize assimilates partitioning while reducing the impact of adverse weather events such as heat stress and frost during the grain-filling stage [START_REF] Camargo | Determining Phenological Patterns Associated with the Onset of Senescence in a Wheat MAGIC Mapping Population[END_REF][START_REF] Chapman | Plant adaptation to climate changeopportunities and priorities in breeding[END_REF][START_REF] Reynolds | Raising yield potential in wheat[END_REF].

Accurate and timely observations of wheat phenology and, particularly, of heading date are, therefore, instrumental for many scientific and technical domains such as wheat ecophysiology, crop breeding, crop management or precision agriculture. Heading date is commonly scored visually in situ by operators that are frequently surveying the crops. It constitutes a labour-intensive task that requires skilled experts. For certain applications like high-throughput phenotyping in crop breeding programs on large genotype panels [START_REF] Araus | Field high-throughput phenotyping: the new crop breeding frontier[END_REF][START_REF] Cabrera-Bosquet | High-throughput Phenotyping and Genomic Selection: The Frontiers of Crop Breeding ConvergeF[END_REF] an accurate visual annotation of the heading date -with experts frequently visiting the field-can be difficult to achieve as the number of microplots increase geometrically. In other applications oriented towards crop management, phenological models may represent a valid alternative to in situ phenology observations [START_REF] Bogard | Predictions of heading date in bread wheat (Triticum aestivum L.) using QTL-based parameters of an ecophysiological model[END_REF][START_REF] White | Simulation-Based Analysis of Effects of and Loci on Flowering in Wheat[END_REF][START_REF] Zheng | Breeding for the future: what are the potential impacts of future frost and heat events on sowing and flowering time requirements for Australian bread wheat ( Triticum aestivium ) varieties?[END_REF], provided that those models are calibrated for each specific genotypes. However, such genotype-specific calibration necessitates extensive phenology observations from field experiments [START_REF] Wallach | How well do crop models predict phenology, with emphasis on the effect of calibration? bioRxiv 708578[END_REF].

The recent development of field sensors and unmanned platforms with imaging capabilities -including aerial (UAVs) and ground vehicles (UGVs)-have opened new avenues to monitor automatically crops in near-real time [START_REF] Baret | Leaf-rolling in maize crops: from leaf scoring to canopy-level measurements for phenotyping[END_REF][START_REF] Comar | A semi-automatic system for high throughput phenotyping wheat cultivars in-field conditions: description and first results[END_REF][START_REF] Jay | Retrieving LAI, chlorophyll and nitrogen contents[END_REF][START_REF] Madec | High-Throughput Phenotyping of Plant Height: Comparing Unmanned Aerial Vehicles and Ground LiDAR Estimates[END_REF][START_REF] White | A Flexible, Low-Cost Cart for Proximal Sensing[END_REF][START_REF] Yang | Unmanned Aerial Vehicle Remote Sensing for Field-Based Crop Phenotyping: Current Status and Perspectives[END_REF]. In parallel, the advances in computational resources and data science achieved in the last years have fostered a significant breakthrough in computer vision and has paved the way to implement advanced algorithms to extract relevant information from high spatial resolution imagery. Deep learning algorithms including convolutional neural network (CNN), have shown excellent performances for object recognition [START_REF] Lecun | Deep learning[END_REF]. These capabilities have favoured their progressive adoption in the fields of agronomy and phenomics [START_REF] Kamilaris | Deep learning in agriculture: A survey[END_REF][START_REF] Singh | Deep Learning for Plant Stress Phenotyping: Trends and Future Perspectives[END_REF]. For instance, they have been successfully used to detect and count individual cereal heads from RGB images [START_REF] Hasan | Detection and analysis of wheat spikes using Convolutional Neural Networks[END_REF][START_REF] Lu | TasselNet: counting maize tassels in the wild via local counts regression network[END_REF][START_REF] Madec | Ear density estimation from high resolution RGB imagery using deep learning technique[END_REF] and LiDAR data [START_REF] Malambo | Automated detection and measurement of individual sorghum panicles using density-based clustering of terrestrial lidar data[END_REF]. Nevertheless, the potential of such algorithms to provide accurate, automatic in situ estimations of crop phenology remains, up to now, underexploited. The use of well-known capabilities of CNNs to detect plant organs on individual images to derive crop phenology from image series needs to be further explored. At the time of writing this article, only very few studies have attempted to do that. The work of [START_REF] Yalcin | Phenology recognition using deep learning: DeepPheno[END_REF] proposing a CNN-based method for the discrimination of phenological stages for several crops, including wheat, seems promising, but no results on the absolute accuracy of the method have been provided. More recently [START_REF] Desai | Automatic estimation of heading date of paddy rice using deep learning[END_REF] have developed a deep learning approach to estimate the heading date on rice, but it was only tested over a small number of situations, preventing from drawing general conclusions about its performance under operational conditions.

The development of operational methods for automatic heading date detection using in situ images would enable to increase enormously the acquisition throughput at a reasonable cost. This would represent an important contribution to the afore-mentioned scientific and technical domains, where frequent phenology observations are critical.

In that context, the objective of this study is to present an automatic method to estimate wheat heading date from daily high-resolution images taken in the field from fixed cameras. The method developed rely on deep learning techniques: daily images are interpreted using a CNN classifier that detects the presence of wheat spikes in the image, and the dynamics of the presence of spikes along the season permits to determine the heading date. A major aspect that differentiates this study against the existing works on this subject is the use of an extensive dataset of observations used to validate the method proposed, which enables to discuss its possible implementation under operational conditions. The dataset comprises 47 field plots sown with several soft and durum wheat cultivars in different regions of France, where fixed cameras were installed, and actual heading dates were annotated by experts. The performance of our CNN-based method to estimate the actual heading date are also concurrently compared with those of a phenology model based on ARCWHEAT [START_REF] Weir | A winter wheat crop simulation model without water or nutrient limitations[END_REF] and calibrated for local conditions. The robustness of our method and its potential operational implementation are discussed with emphasis on possible limitations-e.g. image quality issues, environmental conditions during image acquisition-that may impact the performances. This study was conducted during the years 2017, 2018 and 2019 in different commercial and experimental fields belonging to four contrasted agro-climatic regions around the following cities: Gréoux-les-Bains (43.8° N, 5.9° E) in the south-east of France, Boigneville (48.3° N 2.4° E) in the center of France; Chalons-en-Champagne (49.0° N, 4.4° E) and Saint-Hilaire-en-Woëvre (49.1° N, 5.7° E) in the north-east of France (Figure 1).

Materials and Methods

Study sites

The climate in Greoux-les-Bains is Mediterranean [START_REF] Kottek | World Map of the Köppen-Geiger climate classification updated[END_REF], with a maximum average temperature of 20°C, 690 mm of rainfall. In Boigneville, the climate is temperate and humid, with a maximum average temperature of 15.3⁰C over the year and rainfall of 677 mm (Meteo France). The climate in Chalons-en-Champagne is similar to that of Boigneville, whereas in Saint Hilaire-en-Wovre conditions are slightly colder, especially during winter (Figure 1), and more humid (average precipitation close to 1000 mm/year). Among these four sites, 24 field sensors equipped with RGB cameras (see 2.2) were installed in microplots with a size of 10 x 2 meters belonging to larger experimental fields. The remaining 23 sensors were installed in production plots, with a size similar to a commercial field (around 800 x 200 meters).

Often, the production plots are subdivided in homogeneous units with different cultivars or agro management. In the sites at the north of France, production plots and microplots are sown with winter soft wheat (Triticum aestivum) cultivars Descartes, Oregrain, Fructidor, RGT Sacramento, Matheo and Rebelde. In Greoux-les-Bains, the winter durum wheat (Triticum durum) cultivars RGT Voilur, Anvergur and Toscadou were grown. A summary of the 47 sites considered is given in Table 1. Please refer to Table A in Appendix for a detailed description of the sites, their location and average temperature over the growing season. The field sensors were installed in relatively homogenous areas of the fields which provided daily information over a footprint of about 10m 2 (see 2.2). During the installation, special attention was paid to orientate the camera field-of-view towards the centre of the microplot to prevent possible border effects.

Acquisition of daily images of the canopy with IoTA systems

The field observation systems installed in the 47 sites were developed by Bosch and Hiphen (www.hiphen-plant.com/our-solutions/iot-field-sensor) and are named IoTA (Internet of Things for Agriculture). They consist of a telescopic pole placed vertically and equipped with an RGB camera as well as meteorological sensors (Figure 2a). The RGB camera takes one image (Figure 2b) each day at solar noon, and automatically uploads it to a cloud storage system through a GSM network. The image dimensions are 1024 x 768 pixels and is recorded in PNG format. The camera was set up at a height of approximately 1m above the top of the crop (installed after emergence), with a 45° inclination angle oriented in a compass direction perpendicular to the row. In some of the fields of 2017 and 2018, the length of the pole was adjusted mid-growing campaign to ensure that the camera was always well above the plant canopy. In the 2019 campaign, the height of the pole was fixed at 1.5 m to avoid this midcampaign intervention. It has a field of view (FOV) of 55° x 41° providing a footprint of 10.8 m 2 . Because of the relatively large FOV, the ground resolution is non-uniform throughout the image, particularly in the vertical direction (Figure 2b).

Heading date determined by experts

The heading stage corresponds to the emergence of the developing spike from the flag leaf sheath. [START_REF] Bonnett | The development of the wheat spike[END_REF][START_REF] Zadoks | A decimal code for the growth stages of cereals[END_REF]. In the field, the heading stage is identified according to the definition given by [START_REF] Zadoks | A decimal code for the growth stages of cereals[END_REF]: 50% stems with spikes, at least, half-emerged, corresponding to the phase code 54 of the Zadoks scale. For the 27 sites monitored in 2017 and 2018 (Table 1), an online questionnaire was prepared with IoTA images covering a 9 to 13 consecutive days period approximately centered on the heading date. A panel of 14 experts was asked to identify the reference heading date for all the sites by applying the definition of [START_REF] Zadoks | A decimal code for the growth stages of cereals[END_REF] for heading stage: 50% stems with spikes at least half-emerged. The experts could view the images at their full resolution. Eight of the experts had more than 10 years of experience in working with wheat phenology and only four of them had less than two years of experience. The reference heading date for each of the 27 sites was considered as the average date of all the 14 experts. The standard deviation of the heading dates for each site was also calculated to quantify the variability of the expert replies.

a) b)

In 2019, the actual heading dates were determined in situ by experts on all the 20 available sites (Table 1). They followed the definition by [START_REF] Zadoks | A decimal code for the growth stages of cereals[END_REF] and visited the fields every two or three days.

Note that this reference heading date might be different from the one derived from the 2017-2018 questionnaire since experts were scoring the heading date from images, not from direct observation of the crop in the field.

2.4 Heading date estimation from IoTA images and CNN

Image preparation

To get a more uniform ground pixel size, each image was first cropped into 1024 x 384 pixels by removing the top region of the image where the crop-sensor distance is too large, resulting in blurred objects (Figure 2b). Then 14 overlapping patches of size 256 x 256 pixels were extracted from the cropped image. The overlap between patches was 50% in either vertical or horizontal directions to minimize possible border effects. Working with patches permits to benefit from the full resolution of the IoTA image while avoiding memory issues.

Spikes Labelling

The patches from the 27 sites observed during the 2017 and 2018 growing seasons were labelled into two classes: "spikes present" or "spikes absent". This represents a total of 40,500 patches out of which 17,000 were labelled as "spikes present" and 23,500 as "spikes absent". All the patches belonging to images acquired until five days before the actual heading date determined by experts were automatically assigned to the "spikes absent" class. Similarly, the patches from the images acquired from five days after the actual heading date onwards were assigned to the "spikes present" class. Therefore, only those patches within a window of ±5 days around the actual heading date were visually attributed to their respective classes. Few patches with unclear assignation, such as emerging and sparse spikes were excluded from the training dataset. A few examples belonging to the two classes can be found in Figure 3.

Identifying the presence of wheat spikes with the ResNet50

Wheat spikes were identified in the images using the ResNet50 [START_REF] He | Deep Residual Learning for Image Recognition[END_REF], which obtained the best results in object detection at the ImageNet Large Scale Visual Recognition Challenge 2015 [START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF]. This network has a depth of 50 layers and uses residual blocks with identity mappings. The ResNet50 pre-trained on the ImageNet dataset which is available in the Keras Python deep learning library [START_REF] Chollet | Keras[END_REF] was used.

We replaced the original top layers of the pre-trained ResNet50 by two fully connected layers with dimensions, respectively of 512 and 1 to build a binary classifier (spikes present/absent). The network was re-trained with the labelled image patches by fine-tuning the weights of the entire network with a low learning rate to identify only the high-level features which were relevant to detect spikes and classify the patches as "spikes present" or "spikes absent". This strategy, called transfer learning, performs generally better than training the full network from scratch [START_REF] Lee | Deep-plant: Plant identification with convolutional neural networks[END_REF][START_REF] Tajbakhsh | Convolutional Neural Networks for Medical Image Analysis: Full Training or Fine Tuning?[END_REF]. The • Independent validation on 2019 sites: In this second scheme, a unique set of the 27 sites corresponding to 40500 patches from the 2017 and 2018 sites is used for fine-tuning ResNet50.

The purpose of this scheme is to mimic the operational conditions when CNNs are trained with data from the previous years: the CNN is trained on data acquired in years different from those used for the validation.

In both schemes, 20% of the data were held back for performance testing at the end of each epoch of the ResNet50 re-training. The "binary cross entropy loss" available in the Keras library was used as the loss function. To avoid over-fitting, we reduced the learning rate by a factor of 0.5 when the validation mean absolute error did not improve after three consecutive epochs and stopped the training when the validation mean absolute error did not improve after five consecutive epochs. The pertinence of the ResNet50 network was further evaluated using Gradient-weighted Class Activation Maps (Grad-CAM, developed by Selvaraju et al. , 2016). These maps highlight the regions that contribute to the output score using the gradient values input to the final convolutional layer (shown in Figure 6 and Figure 8).

Heading date estimation

The heading date is determined from the dynamics of the presence of spikes in the images. For each day, , the fraction of patches per image classified by the CNN as 'spikes present' , , is calculated along the growing season (Figure 4). It provides a fair approximation of the Zadoks definition followed by the experts: the proportion of patches with emerged spikes in an image is a reasonable proxy to the proportion of stems with spikes emerged in the image. Then, a three-parameter logistic function is fitted to the time-series of for every site:

Eq. 1

where is the maximum value of , fixed by construction to 1.0 (Figure 4). The maximum growth rate, , and are estimated using the Scipy Python package [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF]. Parameter represents the date when 50% of the patches have spikes.

Heading date from ARCWHEAT model

A version of the phenology module of ARCWHEAT [START_REF] Weir | A winter wheat crop simulation model without water or nutrient limitations[END_REF] adapted to the French local conditions by [START_REF] Gate | Ecophysiologie du ble[END_REF] was used. It is based on cumulated temperature with the effects of vernalization and photoperiod [START_REF] Gouache | Evaluating agronomic adaptation options to increasing heat stress under climate change during wheat grain filling in France[END_REF]. The model was run with actual daily temperatures collected by the Meteo France weather stations that are the closest to the sites in the period 2016-2019.

Cultivar-specific parameters for vernalization and photoperiod corresponding to the nine cultivars observed in this study were adjusted using optimization algorithms based on independent field experiments [START_REF] Thépot | Utilisation d'une population multi-parentale et hautement recombinante de blé tendre pour l'étude de l'architecture génétique de la précocité de floraison[END_REF].

Performance metrics

The accuracy of the CNN classification was evaluated based on the overall accuracy (OA):

Eq. 2 where and are, respectively, the number of patches correctly classified as "spikes present" (true positive) and "spikes absent" (true negative); and ! is the total number of patches in the test dataset.

The root mean squared error (RMSE) was computed to quantify the errors between estimated and observed heading dates:

"#$% & ∑ () * +() , - ./0 1 Eq. 3
where 23 4 is the reference heading date obtained from experts, 23 is the heading date estimated by the indirect method (CNNs or phenology model) and ! is the number of sites used.

Results

Accuracy of heading date estimates from the CNN model

The automatic method proposed estimated the heading date with a RMSE close to 2.0 days (Figure 5), as compared to the reference dates given by the experts. Moreover, the errors when the CNN is trained and validated with images from years 2017 and 2018 using a twofold cross-validation (Figure 5a) are similar to those obtained when validating the method against a completely independent set of images from year 2019 (Figure 5b). Moreover, the coefficient of determination R 2 is very close to 1 especially in 2019, where the variability in the time to heading among the plots monitored is large. In that year, heading was observed about 180 days after sowing for two durum wheat plots (varieties Anvergur and RGT Voilur) at Greoux-les-Bains, more than 20 days earlier than both varieties in 2018. This is explained by higher seasonal temperatures in the 2019 season compared to the 2018 season (seasonal average of 13.3 °C in 2019, against 10.1°C in 2018, see details in the Appendix Table A). standard deviation of the expert annotations for each site in 2017 and 2018. The dotted black lines represent the standard deviation (1.5 days) between experts when identifying the heading dates from photographs.

For two sites at Boigneville monitored in 2018 (Figure 5a), the CNN estimations show discrepancies reaching up to 6 days with the reference dates from the experts. The dynamics of spike appearance for these two sites are shown in Figure 6 along with the GradCAM (Gradient Class Activation Maps) which highlights the regions that influence the CNN prediction. In one of the sites (Figure 6a) the misclassification was due to the poor quality of the images: leaf blades appear bluish because of a poor white balance camera setup that was fixed only after day 203 (Figure 6a). These quality issues introduced substantial artifacts in the time-series that impacted the logistic curve fit. In the second site (Figure 6b), the water droplets on leaves observed on days 195 and 196 after sowing were wrongly identified by the CNN as spikes. Although the CNN seems to slightly over-detect the presence of spikes in this site (Figure 6b) even when no droplets were observed, the errors due to presence of droplets contributed to increase substantially the discrepancies against observed dates. Besides these two specific cases, the issues of CNN misclassification were marginal over the whole dataset. The cross-validation conducted with images from 2017 and 2018 revealed an overall accuracy of 98.45% for classifying individual patches as spikes present/absent. Moreover, the classification errors observed did not exhibit any systematic bias.

Accuracy of heading date prediction by the ARCWHEAT phenology model

The ARCWHEAT crop model, adapted to French conditions had an RMSE of 4.1 days to predict the heading date for the 47 sites monitored from 2017 to 2019 (Figure 7). This represents twice the error yielded by our CNN based method. A similar study conducted by [START_REF] Bogard | Predictions of heading date in bread wheat (Triticum aestivum L.) using QTL-based parameters of an ecophysiological model[END_REF], which evaluated the use of phenology models specifically calibrated for different cultivars to predict heading dates, reported errors comparable to those given by ARCWHEAT in this study. Other works where phenology models are applied at the regional or continental scale with no cultivar-specific calibration show even larger discrepancies reaching up to 20 days [START_REF] Mcmaster | Estimation and evaluation of winter wheat phenology in the central Great Plains[END_REF][START_REF] Ceglar | Improving WOFOST model to simulate winter wheat phenology in Europe: Evaluation and effects on yield[END_REF].

Although the ARCWHEAT model is capable to simulate the variability of the heading date among sites and years (" 5 =0.92), it clearly underperformed when compared to the automatic method based on CNN.

Only 37% of the sites were within the ± 1.5 days interval that represent the variability of the heading dates determined by the experts. Moreover the discrepancies between the model predictions and the

a) b)

observed heading dates were above 10 days for two sites (Figure 7), which is not acceptable for precision agriculture or phenotyping applications.

Figure 7 A comparison between the heading date predictions from the crop phenology model with the reference heading date for all the 47 sites in our study from the 3 growing seasons. The dashed lines around the 1:1 line corresponds to the 1.5 days interval.

To ensure that all the 47 points are visible on the graph despite the overlap, a negligible random noise was added while plotting the points.

Discussion

4.1 Suitability of the automatic CNN-based methodology to estimate wheat heading date in operational applications

Our proposed CNN based method estimated heading dates with about 2 days uncertainties as evaluated over the 47 sites spanning across several regions, years and cultivars. The performances of the method proposed are considered satisfactory since the errors are close to the standard deviation of the expert panel replies provided for years 2017 and 2018. Moreover, this is close to the expected accuracy of an expert visiting the fields every two or three days to annotate crop development stage. The recent study of [START_REF] Desai | Automatic estimation of heading date of paddy rice using deep learning[END_REF] reported estimation errors between one and two days using also a CNN-based approach for paddy rice heading date, but the number of observations and sites were substantially smaller than in our study. Further, among the 20 sites monitored in the 2019 validation dataset, 6 of them were sown with cultivars that were not included in the training set of the CNN, indicating that the method is resistant to possible morphological differences among spikes between cultivars. The robustness of our method is essentially due to three factors: the reliability of the ResNet50 CNN to identify the presence of spikes in the patches; the statistics computed across patches within each image that smooth out individual errors; and the use of a logistic fit to determine the heading date from the daily statistics.

The ResNet50 CNN was used here to classify the patches with spikes. This approach was preferred to directly identifying spikes in the image and counting them [START_REF] Madec | Ear density estimation from high resolution RGB imagery using deep learning technique[END_REF]. Using CNN for image classification increases the efficiency with which the training dataset is generated. Indeed, assigning individual images to the spikes present/absent classes is relatively straightforward: it permitted to generate a training dataset of more than 40,000 patches that contributed to improve the robustness of ResNet50. By contrast, annotating images for spikes identification or spikes counting is time-consuming.

Although CNN-based classification provides much less information than object counting/detection algorithms, the results indicate that this approach is sufficient for heading date estimation. In Figure 8, the first image (top row) spikes are not yet emerged and the Grad-CAM heat map shows low confidence, whereas in the second one (bottom row) the presence of emerged spikes is obvious, and the probability increases up to 0.82. ResNet50 may only detect the presence of a proportion of the spikes present in the patch, but this appears sufficient to correctly classify it. This makes the classification approach more robust for phenology identification as compared to the approaches based on object detection or counting by regression, where the variability of spike size and shape among cultivars or environmental conditions during image acquisition may severely affect the performances [START_REF] Park | Multiresolution Models for Object Detection[END_REF]. The presence of water droplets induced only a moderate bias in heading date estimation in two of the sites (Figure 6b). Images with droplets represent less than 1% of the training dataset cases. That was probably not sufficient to teach the CNN to distinguish between droplets and spikes. These issues can only be solved by increasing the variability by including more images taken under diverse environmental conditions. With 4894 being the prior value of the heading date derived from previous years observations or from phenological models, and ; being the associated maximum error. Anyhow, the use of a logistic function fit makes difficult to estimate the heading date in real time. The correct function fit is only possible a posteriori, i.e. when some images where 1.0 have been already observed, which may only happen 5 to 10 days after the actual heading date. This may limit the use of the method proposed on operational applications of crop management requiring a rapid assessment of the heading date.

Compared to a phenological model calibrated under local conditions, our CNN-based method provides better accuracy while not requiring ancillary information (e.g. sowing dates, variety-specific model parameters) or daily weather data. This is an important advantage, since the accuracy achieved by ARCWHEAT over the 47 sites monitored in this study (RMSE= 4 days) was mainly due to the cultivarspecific calibration of model parameters. However, such a cultivar-specific calibration is time-consuming and requires frequent observations throughout the crop cycle, repeated over several locations and years [START_REF] Cabelguenne | Calibration and validation of EPIC for crop rotations in southern France[END_REF][START_REF] Jégo | Calibration and performance evaluation of soybean and spring wheat cultivars using the STICS crop model in Eastern Canada[END_REF]. Moreover, phenological models largely rely on the quality of meteorological variables, which sometimes are interpolated from weather stations that are far from the sites to be monitored without accounting for possible microclimatic effects [START_REF] Joly | Temperature interpolation based on local information: the example of France[END_REF][START_REF] Monestiez | Spatial interpolation of air temperature using environmental context: Application to a crop model[END_REF]. In any case, mechanistic modelling of plant phenology remains always necessary on prospective studies and to predict cultivar performance under a range of climatic scenarios. It also constitutes an attractive alternative for heading date estimation for some applications where the acquisition of high-throughput canopy images is not feasible. Phenological models can also largely benefit from the method developed in this paper: the use of IoTA field systems and deep learning approaches would substantially reduce the cost of calibration experiments and would permit also to increase the environmental variability of the field trials giving access to frequent observations from remote locations.

Our CNN based method constitutes a robust and cost-efficient approach for heading date estimation for operational applications when daily images of the canopy are available, as it is the case for some highthroughput phenotyping platforms. In those platforms, vectors such as unmanned ground vehicles and hand-pushed carts are often used to frequently monitor the plant development and characterize biophysical traits of the different cultivars using optical images [START_REF] Deery | Proximal Remote Sensing Buggies and Potential Applications for Field-Based Phenotyping[END_REF][START_REF] Mueller-Sim | The Robotanist: A ground-based agricultural robot for high-throughput crop phenotyping[END_REF][START_REF] White | A Flexible, Low-Cost Cart for Proximal Sensing[END_REF]. In such applications, the proposed automatic method showcased in this study could be directly integrated into other data processing pipelines at almost no cost to estimate heading date from the existing RGB images.

In agronomic or breeding applications where the phenology on distant fields need to be monitored -e.g.

regional or national networks of agricultural fields-the whole system presented, including the fixed camera and the CNN-based method to estimate heading dates, may present important advantages in terms of costs compared to in situ visual annotation. The initial investment in each camera (including the pole, batteries and the hardware for data transmission) raises up to, approximately, 650 €; and the yearly costs of system maintenance and real-time data transmission by GSM is about 150 € per camera.

The cost-efficiency -and the benefits in terms of environmental impact-of a system based on network of cameras compared to expert visits has to be determined case by case, and will largely depend on the distance between the fields to be monitored: the larger the network is, the more efficient remote observations are compared to field visits.

Subjectivity of visual annotations of the heading date from experts using RGB images

The visual determination of heading date is a task that includes some degree of subjectivity. For the 27 sites monitored in 2017 and 2018, the panel of 14 experts provided different estimates of the heading dates after inspecting visually the IoTA images. In most of the cases, the panel members selected 5 to 7 different dates per site out of those proposed in the questionnaire, although they were asked to follow the same definition of the heading date [START_REF] Zadoks | A decimal code for the growth stages of cereals[END_REF]. The distribution of the deviations between individual replies and the average date for each site was roughly Gaussian with a standard deviation of 1.5 days (Figure 9), very close to the RMSE of our CNN based method proposed in this paper.

The subjectivity when determining visually heading date from images is obviously higher compared to in situ scoring. Issues regarding the image quality exacerbated the discrepancies among experts. According to the feedback provided by the panel, image saturation and blur due to a suboptimal camera setup made difficult to see the emergence of spikes in some of the images. Further, spikes were harder to identify in images taken under direct illumination conditions due to low image contrast, as well as for cultivars with awns.

The variability of expert replies evaluated in the 27 sites constitutes a good benchmark for our CNN based method, as image quality issues may also affect the identification of spikes by ResNet50 when classifying patches. The similarity between the 2 days RMSE associated to our CNN based method with the 1.5 days confidence interval of the expert date demonstrates that the performances of our CNN based method can be considered comparable to the expert reply when observing the same images.

Conclusion

In this study, we propose a CNN based method to estimate the heading date from daily images acquired over wheat crops using an RGB camera fixed in the field. Images are processed per patch based on a binary CNN-classifier to construct the dynamics of spikes appearance. Our method is easy to implement since the labelling of patches is not time-consuming as compared to individual object annotation required for other CNN models used for object identification or counting. The reliability of the CNNclassifier to identify the presence of emerged spikes was marginally affected by the illumination conditions and cultivar diversity, since the training dataset included images acquired under diverse environmental conditions. Our method achieved satisfactory performances with RMSE≈2.0 days, which is close to the uncertainties of expert annotations, and substantially better than phenological models specifically calibrated for the cultivars monitored.

The proven robustness of the proposed method suggests a strong potential for cost-efficient operational applications in the field of phenotyping and agronomic applications. However, our method is limited by the image footprint close to 5 m² if only the bottom half of the image is used, as done in our study.

However, the good consistency with expert observations taken over a larger sampling area demonstrates that our restricted sampling was sufficient. Nevertheless, the representativeness of such small footprint estimations to characterize phenology over large and heterogenous fields remains an open question for future works.

This method has been developed using a camera looking to the crop in a fixed position, but a similar approach could be transposed to time series of images from other vectors used in phenotyping experiments, such as unmanned ground and aerial vehicles, providing that the revisit time and resolution are sufficient. Further, the method could be adapted to identify other crop development stages associated with the identification of certain organs, such as the appearance of anthers for wheat to date flowering, or the appearance of tassels for the male flowering in maize.

sugar beet crops from multi-angular optical remote sensing:

Figure 1 .

 1 Figure 1. Location of the 4 regions where field sensors installed during the 2017, 2018 and 2019 growing seasons. The graph indicates the average monthly temperatures in the period 1986-2018 observed in 4 near weather stations from INRA and Meteo france: St. Gilles (Gréoux-les-Bains), Gif-sur-Yvette (Boigneville, only since 1999), Erneville (St-Hilare-en-Wovre), and Fagnieres (Chalons-en-Champagne).

Figure 2

 2 Figure 2 a) The IoTA systems installed in a field. The top part inclined at 45° hosts the RGB camera. The cylinder attached to the vertical pole includes sensors to monitor the local air temperature and moisture. b) An example of a raw image.

Figure 3

 3 Figure 3 Samples extracted from patches belonging to the 9 wheat varieties monitored in our study. The patches were classified as 'spikes present' (top row) and 'spikes absent' (bottom row).

  training dataset comprises the 27 sites available in 2017 and 2018. First-order data augmentation was applied to the patches: translation, rotation, zoom, flip and changes in brightness levels. This improves the generalization capacity of the neural network and increases the size of the training set at marginal cost. Two different training and validation schemes were followed:• Twofold-cross validation on 2017 and 2018 sites: the 27 sites available in 2017 and 2018 were divided randomly into folds of 13 and 14 sites with approximately 19500 and 21000 labelled patches in each fold. ResNet50 was re-trained on each of the two folds independently and validated with the other one.

Figure 4

 4 Figure 4 Typical dynamics of the fraction of patches with spikes ( shown as green dots). The logistic curve fitted to the dynamics of is shown using dashed blue line. The heading date estimated from adjusted parameter is represented by the vertical blue line.

Figure 5

 5 Figure 5 Comparison of reference heading dates and those estimated by our CNN-based method. a) 2017 and 2018, using a twofold cross-validation for years 2017 and 2018; b) independent validation over year 2019. Horizontal error bars represent the

Figure 6

 6 Figure 6 Dynamics of for two sites in Boigneville in 2018 where large errors on heading date estimations by the CNN method were observed due to non-optimal image quality (a); image misclassification due to the presence of water droplets (b).

Figure 8

 8 Figure 8 Patches extracted from two images of the same IoTA system and acquired on the reference heading date (top) and one day after (bottom) with their corresponding class activation maps. The confidence score provided by the CNN towards the presence of spikes is shown along with the Grad-CAM heat map which shows the areas which influence the CNN output

Figure 9

 9 Figure 9 Distribution of the difference in days between the individual heading date identified by the experts and the mean heading date value.

Table 1

 1 Description of the distribution of the 47 sites across years and regions. The number of sites available each year per region and year is indicated along with the corresponding number of cultivars.

			2017		2018		2019
	Locations	Sites	Cultivars	Sites	Cultivars	Sites	Cultivars
	Gréoux-les-Bains	-	-	7	3	2	2
	Boigneville	8	1	12	3	12	3
	Chalons en Champagne	-	-	-	-	3	1

  Comparison of vegetation indices and

	PROSAIL	inversion	for	field	phenotyping.	F.	Crop.	Res.	210,	33-46.

Acknowledgements

The authors would like to thank the Association Nationale de Recherche Technologique (ANRT) for their grant CIFRE. The authors acknowledge Guy Deshayes (Arvalis) for his kind support with consolidating the reference data from the experts in the field.


Chloé Roullier-Gall
email: chloe.roullier-gall@u-bourgogne.fr

Julie Signoret

Christian Coelho

Daniel Hemmler

Mathieu Kajdan

Marianna Lucio

Bernhard Schäfer

Régis D Gougeon

Philippe Schmitt-Kopplin
email: schmitt-kopplin@helmholtz-muenchen.de

Philippe Schmitt- Kopplin

Chloé Roullier

Influence of regionality and maturation time on the chemical fingerprint of whisky

Keywords: Whisky, metabolomics, chemometrics, authentication, geographical origin and maturation time

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Across Scotland there are up to 100 operating distilleries producing whisky distributed in five main regions commonly recognized as Lowlands, Highlands, Speyside, Islands (including whisky from Orkney, Skye, Mull, Jura, Arran and the Campbeltown peninsula) and Islay [START_REF] Barnard | The Whisky Distilleries of the United Kingdom[END_REF].

The composition of whisky is complex and mainly impacted by the water and the cereals used, the fermentation, the distillation, the maturation and the blending process [START_REF] Garcia | Whisky analysis by electrospray ionization-Fourier transform mass spectrometry[END_REF][START_REF] Piggott | Science and technology of whiskies[END_REF][START_REF] Reid | Assessment of Scotch whisky quality by pyrolysis-mass spectrometry and the subsequent correlation of quality with the oak wood cask[END_REF].

Traditionally, whiskies are matured in oak wood barrels for a minimum of three years before consumption, but many whiskies are matured for 12 or more years(JOHN R. [START_REF] Piggott | Effects on Scotch whisky composition and flavour of maturation in oak casks with varying histories[END_REF]. As for wine, sensory properties of whiskies are known to improve during the period of maturation in the oak barrel, where the composition changes through a complex array of chemical reactions [START_REF] Macnamara | Advances in the Ageing Chemistry of Distilled Spirits Matured in Oak Barrels[END_REF][START_REF] Macnamara | Flavour components of whiskey. III. Ageing changes in the low-volatility fraction[END_REF]. Previous studies have attempted to describe the chemical changes during maturation, with particular emphasis on the color, the aroma, the volume, the strength and the flavor compounds (Lee, Paterson, Piggott, & Richardson, 2001a;[START_REF] Liebmann | Changes in whisky while maturing[END_REF][START_REF] Macnamara | Advances in the Ageing Chemistry of Distilled Spirits Matured in Oak Barrels[END_REF][START_REF] Macnamara | Flavour components of whiskey. III. Ageing changes in the low-volatility fraction[END_REF].

Due to its large commercialization and relatively high prices, whisky counterfeiting and adulteration is quite common worldwide [START_REF] Garcia | Whisky analysis by electrospray ionization-Fourier transform mass spectrometry[END_REF][START_REF] Wiśniewska | Chemical composition analysis and authentication of whisky[END_REF]. In case of fraudulent Whisky, maturation time and history of the oak casks can be mislabeled [START_REF] Stupak | A novel approach to assess the quality and authenticity of Scotch Whisky based on gas chromatography coupled to high resolution mass spectrometry[END_REF]. Authenticity is one of the major concerns for the distillers, dealers and consumers of whiskies around the world(R. I. [START_REF] Aylott | Analytical strategies to confirm the generic authenticity of Scotch whisky[END_REF][START_REF] Aylott | Analytical strategies to confirm Scotch whisky authenticity[END_REF][START_REF] Heller | A rapid sample screening method for authenticity control of whiskey using capillary electrophoresis with online preconcentration[END_REF]; M. MacKenzie & I. [START_REF] Mackenzie | Analytical strategies to confirm Scotch whisky authenticity. Part II: Mobile brand authentication[END_REF][START_REF] Møller | Electrospray ionization mass spectrometry fingerprinting of whisky: immediate proof of origin and authenticity[END_REF][START_REF] Parker | Investigation into the use of carbon isotope ratios (13C/12C) of Scotch whisky congeners to establish brand authenticity using gas chromatography-combustion-isotope ratio mass spectrometry[END_REF][START_REF] Stupak | A novel approach to assess the quality and authenticity of Scotch Whisky based on gas chromatography coupled to high resolution mass spectrometry[END_REF]. Both, maturation period and history of the casks in which maturation occurred are important for the final composition [START_REF] Kew | Chemical Diversity and Complexity of Scotch Whisky as Revealed by High-Resolution Mass Spectrometry[END_REF][START_REF] Roullier-Gall | Usage of FT-ICR-MS Metabolomics for Characterizing the Chemical Signatures of Barrel-Aged Whisky[END_REF]. Recently, Whisky authentication has been the subject of several studies employing various analytical techniques, such as: GC and LC with different detectors (FID, AED, UV-Vis), electronic nose, atomic absorption spectroscopy and mass spectrometry [START_REF] Wiśniewska | Chemical composition analysis and authentication of whisky[END_REF]. The aroma and flavor profiles are classically used to distinguish whiskies by distillery or even to reflect their places of origin [START_REF] Jack | Modelling the sensory characteristics of Scotch whisky using neural networks-a novel tool for generic protection[END_REF][START_REF] Jackson | Malt Whisky Companion[END_REF][START_REF] Lee | Sensory discrimination of blended Scotch whiskies of different product categories[END_REF]. Gas chromatography coupled to tandem mass spectrometry (GC-Q-ToF) was employed to distinguish malt whiskies according to the type of cask in which they were matured (bourbon versus bourbon and wine) [START_REF] Stupak | A novel approach to assess the quality and authenticity of Scotch Whisky based on gas chromatography coupled to high resolution mass spectrometry[END_REF]. Infrared spectroscopy with statistical analysis allowed to distinguish Scottish, Irish, and American whiskies and 2 and 3 years old beverages from 6 and 12 years old whiskies [START_REF] Sujka | The application of FT-IR spectroscopy in discrimination of differently originated and aged whisky[END_REF]. Recently, two different, hypothesis-free, sensor arrays based upon three fluorophores successfully discriminated whisky samples with respect to origin (American, Irish and Scotch Whisky) and taste (rich vs. light) [START_REF] Han | A Hypothesis-Free Sensor Array Discriminates Whiskies for Brand, Age, and Taste[END_REF]. Mass spectrometry is more and more used for quality control and proof of authenticity of whisky samples. Direct infusion electrospray ionization mass spectrometry (ESI-MS) was used in order to provide a direct, rapid and sensitive method for the characterization of distilleries and authenticity of whisky samples [START_REF] Møller | Electrospray ionization mass spectrometry fingerprinting of whisky: immediate proof of origin and authenticity[END_REF]. Garcia et al., highlighted the use of Fourier-transform ion cyclotron resonance mass spectrometry (FT-ICR-MS) as an approach to screen for ageing and counterfeiting [START_REF] Garcia | Whisky analysis by electrospray ionization-Fourier transform mass spectrometry[END_REF]. FT-ICR-MS was used by [START_REF] Kew | Chemical Diversity and Complexity of Scotch Whisky as Revealed by High-Resolution Mass Spectrometry[END_REF] for the discrimination of whiskies according to their blending process and type of casks used in maturation [START_REF] Kew | Chemical Diversity and Complexity of Scotch Whisky as Revealed by High-Resolution Mass Spectrometry[END_REF]. More recently, non-targeted FT-ICR-MS combined with LC-MS/MS showed the impact of the wood and in particular the history of the barrel on the distillate composition during ageing and revealed the importance of the initial composition of the distillate and the distillery process [START_REF] Roullier-Gall | Usage of FT-ICR-MS Metabolomics for Characterizing the Chemical Signatures of Barrel-Aged Whisky[END_REF].

The objective of this work was to go beyond our previous work and implement strategies by using statistical modeling to discriminate whiskies according to their geographical origin and authenticate the maturation time in the barrel based on the non-volatile profiles. Here, we show that the combination of FT-ICR-MS and statistical analysis of whiskies from 3 to 43 years of cask maturation allowed the authentication of the geographical origin and maturation time. The combination of metabolomics and chemometrics was able to reveal chemical markers associated with the ageing, regardless of the origin or the production process.

Materials and Methods

Whisky samples

106 whisky samples from 32 different distilleries in Scotland and from 3 to 43 years of maturation were analyzed (Supplemental table 1). All samples were collected directly from the bottle and stored in 10 mL amber vials at room temperature until analysis.

Direct infusion FT-ICR-MS

Ultrahigh-resolution FT-ICR mass spectra were acquired with a 12 T Bruker SolariX mass spectrometer (Bruker Daltonics, Bremen, Germany) equipped with an APOLLO II electrospray source in negative ionization mode. For MS analysis, whiskies were diluted 5:100 (v/v) with methanol (LC-MS grade, Fluka, Germany). The diluted samples were infused into the electrospray ion source with a flow rate of 120 μL h -1 . Operating parameters of direct infusion FT-ICR-MS were carried out according to [START_REF] Roullier-Gall | Usage of FT-ICR-MS Metabolomics for Characterizing the Chemical Signatures of Barrel-Aged Whisky[END_REF] previous study. Settings for the ion source were: drying gas temperature 180 °C, drying gas flow 4.0 L min -1 , capillary voltage 3,600 V. Spectra were first externally calibrated by ion clusters of arginine (10 mg mL -1 in methanol). Internal calibration of each spectrum was conducted with a reference list including selected whisky markers and ubiquitous fatty acids.

The spectra were acquired with a time-domain of 4 megawords and 400 scans were accumulated within a mass range of m/z 92 to 1000 corresponding to a run time of 18 min. A resolving power of 400,000 at m/z 300 was achieved. Quality control (QC) samples were prepared by pooling equal amounts of all samples. QC samples were analyzed at the beginning and after every tenth sample to monitor the reproducibility of the measurements [START_REF] Roullier-Gall | Usage of FT-ICR-MS Metabolomics for Characterizing the Chemical Signatures of Barrel-Aged Whisky[END_REF].

Processing of FT-ICR-MS data

Raw spectra were post-processed by Compass DataAnalysis 4.2 (Bruker Daltonics, Bremen, Germany) and peaks with a signal-to-noise ratio (S/N) of at least 6 were exported to mass lists. All exported m/z features were aligned into a matrix containing averaged m/z values (peak alignment window width: ±1 ppm) and corresponding peak intensities of all analyzed samples. Molecular formulae were assigned to the exact m/z values by mass difference network analysis using an in-house developed software tool. In total, 5979 detected features could be assigned to distinct and unique molecular formulae. More than 90% of all assignments were found within an error range lower than 0.2 ppm. All further calculations and filtering were done in Perseus 1.5.1.6 (Max Planck Institute of Biochemistry, Germany) and R Statistical Language (version 3.1.1).

Repeatability of FT-ICR-MS whisky measurements

The repeatability of the measurements was evaluated on QCs samples injected in triplicate at the beginning and at the end of the sequence analysis. The repeatability over time of FT-ICR-MS analyses during the entire runtime was also evaluated , by monitoring the precision of the recorded peak intensities from ten analyzed QC samples and manual comparison of the raw mass spectra [START_REF] Roullier-Gall | Usage of FT-ICR-MS Metabolomics for Characterizing the Chemical Signatures of Barrel-Aged Whisky[END_REF]. In the QC samples, 92.8% of all detected monoisotopic signals showed a relative standard deviation lower than 20% [START_REF] Roullier-Gall | Usage of FT-ICR-MS Metabolomics for Characterizing the Chemical Signatures of Barrel-Aged Whisky[END_REF].

Excitation Emission Matrix Fluorescence

Whisky fluorescence was acquired on a spectrofluorometer (Horiba Aqualog). Excitation ranged from 225 to 600 nm and emission was recorded from 230 to 600 nm (C. [START_REF] Coelho | Fluorescence Fingerprinting of Bottled White Wines Can Reveal Memories Related to Sulfur Dioxide Treatments of the Must[END_REF]. The spectrofluorometer used a conventional right-angle optical setup. Excitation Emission Matrices (EEM) were corrected daily for their Rayleigh and Raman scattering, as well as inner filtering effects by using the functions provided within the FluorEssence software of the Horiba Aqualog system. The fluorescence intensity was normalized to a Starna 1 ppm quinine sulfate reference cell, measured every day.

UPLC targeted analysis of polyphenols

An Acquity Waters ultra performance liquid chromatography (UPLC) with a diode array fluorescence detector (DAD) was used to separate and quantify individual polyphenolic compounds in whisky samples. The column was a BEH C18, 1.7 µm, 2.1 mm × 150 mm.

Column temperature was kept constant at 30 °C, and samples were held at 8 °C. An elution was applied, starting isocratic from 100% A (ultrapure water, 0.1% formic acid) from 0 to 6 min and then the gradient was increased linearly over 56 min to 100% B (methanol, 0,1% formic acid), where it was held again until reaching 60 min, with a rate of 0.3 mL•min -1 .

Injection volume was 5 μL. Detection limits (LODs), quantification limits (LOQs) and intraand inter-day variation of the assay were determined by injecting a series of dilute solutions with known concentrations.

Total phenolic content (TPC) of whiskies

Total phenolic content (TPC) in whisky was determined using the Folin-Ciocalteu assay: 0.1 mL of undiluted whisky, 0.5 mL of Folin-Ciocalteu reagent, 2 mL of a 20 (w/v) sodium carbonate solution and 7.4 mL of ultrapure water were mixed and reacted for 30 min at room temperature. Then, the absorbance at 750 nm was measured and calculated as mg.L -1 of gallic acid equivalent using a calibration curve between 0 and 500 mg.L -1 of gallic acid. All samples were analyzed in triplicates.

Data analysis

The statistical analysis was performed with Perseus, Simca 12.0 (Umetrics, Umeå, Sweden) and R (Version 1.0.136).

Van Krevelen and elemental diagrams were used to visualize FT-ICR MS data [START_REF] Gonsior | Photochemically Induced Changes in Dissolved Organic Matter Identified by Ultrahigh Resolution Fourier Transform Ion Cyclotron Resonance Mass Spectrometry[END_REF][START_REF] Gougeon | Authentication Approach of the Chemodiversity of Grape and Wine by FTICR-MS[END_REF][START_REF] Hertkorn | Natural Organic Matter and the Event Horizon of Mass Spectrometry[END_REF]. Multivariate statistical methods (.

Principal component analysis and hierarchical cluster analysis) were used to explore similarities and hidden patterns among samples [START_REF] Barker | Partial least squares for discrimination[END_REF][START_REF] Granato | Use of principal component analysis (PCA) and hierarchical cluster analysis (HCA) for multivariate association between bioactive compounds and functional properties in foods: A critical perspective[END_REF]. As within-groups variability dominates the among-groups variability for the geographical origin and maturation time impact, partial least square discriminant analysis were used [START_REF] Barker | Partial least squares for discrimination[END_REF]. We set up diverse classification models in order to reveal the different geographical origin and maturation time for barrel specific metabolites. In order to improve the efficiency of such classification and the possible presence of overfitting and noise, we preprocessed the entire dataset applying the ReliefF algorithm [START_REF] Witten | Data Mining: Practical Machine Learning Tools and Techniques[END_REF]. The algorithm identified a subset of variables that was used for the next classification models maximizing the model performance in term of accuracy. The features' selection was based on the highest rank value attributed to each variable by the algorithm. Consequently, with the reduced datasets we built several partial least square discriminant analysis (PLS-DA) models. The goodness of the fit was evaluated by the coefficient of determination (R²X), the proportion of the variance of the response that is explained by the model (R²Y), the predictive ability (Q²Y) and validation was carried out using permutation tests. For the ageing time model, whisky samples were randomly divided into two groups: 101 samples selected to build the statistical model and the 5 remaining samples employed as external validation. A Bland-Altman plot was used to analyze the agreement between the known ageing time and the model prediction.

Results and discussion

Geographical origin impact of Scotch whisky

A visual comparison of FT-ICR-MS spectra of Scotch whiskies from Highland, Lowland, Speyside and Islay distilleries and from 3 to 43 years of cask maturation is presented in and show comparable intensities in all four whisky samples (Figure 1). The Venn and van Krevelen diagrams (Figure 1) show that 2200 formulas (out of 5979) were found in all samples of the four groups of whisky samples confirming a great consistency in the chemical composition regardless of the origin [START_REF] Kew | Chemical Diversity and Complexity of Scotch Whisky as Revealed by High-Resolution Mass Spectrometry[END_REF].

Despite the high similarities between whiskies (36.8% of formulas were found in all samples of the four groups of whisky samples) as shown in the Venn diagrams and the direct comparison of spectra (Figure 1), 533 (8%), 610 (10%), 401 (6.7%) and 464 (7.7%) formulas were uniquely found in at least 50% of Highland, Islay, Speyside and Lowland whisky samples, respectively. In addition to the unique formulas, differences in signal intensities between whiskies can be found which could result from particular distilleries, maturation times, geographical features or environmental impacts. As example, the signal at m/z 261.06158 corresponding to the molecular formula [C10H13O8] -is less intense in the whisky from Speyside (2•10 7 ) compared to the other three whiskies (5-6•10 7 ) (Figure 1).

The impact of the type of whisky (blend or malt), the maturation wood type and the wood history on the Whisky composition has recently been studied [START_REF] Kew | Chemical Diversity and Complexity of Scotch Whisky as Revealed by High-Resolution Mass Spectrometry[END_REF][START_REF] Roullier-Gall | Usage of FT-ICR-MS Metabolomics for Characterizing the Chemical Signatures of Barrel-Aged Whisky[END_REF]. However, barley, peat and water origin, malting, brewing, fermentation and distillation process, nature of the cellars, climates as well as environmental factors could also impact the complex composition of whisky. Figure 2 shows the first two principal components obtained from principal component analysis computed on the 106 analyzed whisky samples from Scotland. Whiskies showed a rough separation by their maturation time in the barrel on principal component 2. By comparison, the impact of the geographical origin -Highlands, Lowlands, Speyside and Islay -on the whisky metabolome appears less significant in the PCA compared to the maturation time (Figure 2). The samples from the Highlands and Lowlands appeared in a more discrete cluster than the Speyside and Islay samples, which are more distributed all around the plot. Interestingly, the four samples in the top left corner of the PCA were from the same distillery (but with different maturation times) and the eight samples (in red) in the top right corner were from different distilleries but all in Islay and all with six years of barrel maturation (Figure 2). The great diversity in the composition of the whisky samples demonstrates the large number of environmental parameters (including production, geographical origin, maturation time…), which eventually impacts the final product.

A group separation is not easily observable in PCA models because of the high number of sample characteristics (maturation from 6 to 43 years, 4 geographical origins and 32 manufacturers) and the high number of variables considered (5979 molecular formulas). A PLS regression model was used to better distinguish the geographical origins of the whisky samples. Figure 3 shows results from a PLS regression computed from FT-ICR-MS data of 106 whisky samples from the four areas in Scotland, 32 different distilleries and from 3 to 43 years of maturation. The first two components allowed to separate whiskies into four distinct groups according to the geographical origin (Highlands, Lowlands, Speyside and Islay) and independent of the distillery or the maturation time (Figure 3). Whiskies from Lowlands and Islay distilleries appear to be chemically more different from the samples from Highlands and Speyside distilleries (separation on the 1 st component). The second component allowed to distinguish whiskies from the Highlands and Lowlands from those produced in Speyside and Islay distilleries. Based on the PLS, specific masses for each location were extracted (VIP masses) and projected into van Krevelen diagrams (Figure 3). For example, samples from Speyside clearly showed an enrichment of polyphenols or polyphenol-type compounds (Figure 3). By comparison, Highland samples seem to be characterized by higher alcohols with a profile reminiscent to the distillate [START_REF] Roullier-Gall | Usage of FT-ICR-MS Metabolomics for Characterizing the Chemical Signatures of Barrel-Aged Whisky[END_REF]. Lowland and Islay appear to be characterized by carbohydrates. The schematic classification of single malts according to their origin could be based, in part, on a geological reality, of which water would be the main factor. Meir-Augenstein showed that the source water used for and during whisky production correlates with the authenticity of whisky(Meier-Augenstein, Kemp, & Hardie, 2012). Water is involved in several stages of processing: added to ground barley to produce the mash, when diluting the distilled alcohol to cask strength for maturation and when reducing cask strength whisky for bottling. Meir-Augenstein's study suggests that whisky keeps a signature of the geographic provenance of the used water [START_REF] Meier-Augenstein | Detection of counterfeit scotch whisky by 2H and 18O stable isotope analysis[END_REF]. However, from all elements necessary for whisky production, water is the one whose impact on the final composition is probably the most difficult to evaluate [START_REF] Witten | Data Mining: Practical Machine Learning Tools and Techniques[END_REF].

Another source of the variability in the chemical composition is the distillery brand.

Compositional specificities from manufacturers can contribute to the geographical region model. Studying samples from Islay and from three distilleries including various maturation times highlights the ability of FT-ICR-MS to differentiate whiskies according to the distilleries and independently of the maturation time. Non-supervised statistics allowed to visualize the impact of the process, composition and factory on the final whisky independently of the maturation (supplementary Figure 1). According to the hierarchical cluster analysis (supplementary Figure 1), whiskies from the manufacturers W2 and W1 manufacturers seem close and distinct from W3 although W2 and W3 are geographically very close to one another (less than 3 km). According to the van Krevelen diagrams in supplementary Figure S1, the whisky from distillery W1 (including whiskies from 12 to 18 years in barrel) seems to be characterized by a higher concentration in polyphenolic compounds whereas whiskies from distillery W2 (12 to 21 years in barrel) are mainly characterized by carbohydrate-type compounds as well as sulfur and nitrogen containing compounds. The whisky from distillery W3 (6 to 12 years in barrel) seems to be characterized by higher alcohols and shows a profile reminiscent to distillate [START_REF] Roullier-Gall | Usage of FT-ICR-MS Metabolomics for Characterizing the Chemical Signatures of Barrel-Aged Whisky[END_REF].

The samples differing in brand and geographical origin can be properly grouped and characterized by FT-ICR-MS and chemometric analysis. These combined techniques enable the fast evaluation and authentication of Scotch whiskies independently of the maturation time or the process used. They allowed an excellent determination of whisky composition as well as distinguishing between different brands of whiskies produced in areas only 20 km apart from one another.

Impact of maturation time on the chemical composition of Scotch whiskies

Based on the huge impact of maturation in wood barrels on the chemical diversity of whisky, as shown in our previous study [START_REF] Roullier-Gall | Usage of FT-ICR-MS Metabolomics for Characterizing the Chemical Signatures of Barrel-Aged Whisky[END_REF], we here focused on the impact of the maturation time on the final composition. We first studied three samples from one distillery (W3) in Islay after 6, 8 and 10 years of cask maturation (Supplementary Figure 2).

The samples revealed 376 mass peaks whose intensity significantly increased with maturation time, and 81 mass peaks with decreasing intensity (Supplementary Figure 2a-b). The higher number of features which showed an increase in peak intensity confirms that chemicals transfer from wood to the distillate as a function of ageing.

In order to extract age-related metabolite signatures, a PLS discriminant analysis was designed using 7 groups with different maturation times: whiskies that had been aged in casks for 6 years, 8 years, 10 years, 12 years, 14 years, 16 years and 18 years (Supplementary Table 3). We selected 76 from the 106 analyzed whisky samples which ranged from 6 to 18 years.

A PLS-DA model was built based on 71 of these whisky samples (Figure 4c). Five of the 76 samples were used for model validation (Figure 4e). The separation of whiskies according to the ageing time (Figure 4c) showed a good predictive power with Q² of 0.49 and R²Y of 0.77.

The scatter plot of the model group versus the predicted group for each sample confirmed the robustness of the prediction model (Figure 4 D), with a prediction average close to the maturation time and an excellent standard deviation (between 0.69 to 1.57 years, Supplementary Table 3). However, it was not possible to really differentiate groups of 12 and 14 years of maturation. These samples were predicted as 12.3 years for both, the group of 12 years and 14 years old whiskies (maximum standard deviation was 1 year).

Finally, the five whisky samples, which were not used for the model construction (maturation of 6, 10, 12, 12 and 16 years in barrel) were used to validate the model. The good predictions, achieved for these five whisky samples, confirmed the predictive ability of the maturation time model (Figure 4 C andE). After model validation, 256 discriminating and characteristic masses (VIP variables higher than 1), whose peak intensities increased with maturation, were extracted (Figure 4 F and H). As example, compounds from wine such as syringic acid and caftaric acid and from wood such as lyoniresinol, patuletin and digallic acid were found as increasing with maturation time (Figure 4 and supplementary Figure 2). In contrast, 213 masses with decreasing intensity led to 54 annotations from databases including stearyl acetate and syringaledehyde (Figure 4).

ANOVA tests were applied to extract subsets of features that discriminate between whiskies with a maturation time of 6 and 8 years (supplementary figure 4A), 8 and 10 years (supplementary figure 4B), 10 and 12 years (supplementary figure 4C), 12 and 14 years (supplementary figure 4D), 14 and 16 years (supplementary figure 4E), and 16 and 18 years (supplementary figure 4F), respectively. Statistically extracted masses can be used to highlight the impact of maturation time on the whisky composition between the different ageing groups. According to the statistics, 27 molecular formulas were significantly enriched in the youngest whiskies (6 years of maturation) while 36 formulas were specific to the 66 whiskies older than 6 years of maturation. In the same way, 74 molecular formulas were statistically identified as markers to differentiate between whiskies younger than 8 years of maturation (6 and 8 years of maturation) and 132 formulas were identified as markers for whiskies older than 8 years of maturation. The number of significant masses statistically extracted increased until the distinction of whiskies between 10 (121 molecular formulas) and 12 (1182 molecular formulas) years of casks maturation. Then the number of significant masses starts to decrease.

Interestingly, the higher differences were found between whiskies after 10 and 12 years of maturation directly followed by whiskies after 14 and 16 years of maturation (supplementary Figure 4). The similarities between samples characterized by 12 and 14 years of maturation, which probably can be explained by the low ability of the model prediction to distinguish those two groups (Figure 4). The difference between young (less than 12 years of maturation, group 1) and old whiskies (more than 12 years of maturation, group 2) was further studied by Excitation Emission Matrix (EEM) fluorescence spectroscopy and analysis of polyphenols (supplementary figure 5 and6). As illustrated in Supplementary figure 6A, PARAFAC results enabled to easily discriminate whiskies from both groups on the first component.

Whiskies of group 1 (younger than 12 years of maturation) were driven by wood barrel ageing processes, conferring a relative homogeneity to this group. Whiskies from the group 2 (older than 12 years of maturation) showed a higher within-group-variation (confirming results from FT-ICR-MS), which may be because of different cask finishing practices. This may include re-casking for a secondary maturation into other wood casks that were used for ageing of wines, spirits or even beers (Scotch Whisky Association, n.d.). Such practice tends to diversify Scottish whisky profiles by masking the original distillate. The fluorescence positions of the PARAFAC components suggest differences in polyphenolic compounds (Figure 6B) which could be confirmed in a TPC assay and quantification of syringic acid, gallic acid and scopoletin by UPLC (Figure 6C). Scopoletin was also identified as wood marker increasing in intensity with maturation time. No significant differentiations were found for tyrosol between our two groups, surely due to a moderated retention of this compound to wood surfaces (Christian [START_REF] Coelho | Chemical Transfers Occurring Through Oenococcus oeni Biofilm in Different Enological Conditions[END_REF].

Out of the 1182 molecular formulas found significantly increased for whiskies with at least 12 years of maturation, only 222 annotations were found in databases: 167 were found in our home-built plant database, 154 in Metlin, 130 in KEGG and 25 in HMDB (Supplementary Figure 7 and supplementary Table 4). Interestingly, 467 of the 1182 markers could be previously identified as wood specific markers (Roullier-Gall et al., 2018) (Supplementary Figure 8). Moreover, 675 of the 1182 features could be detected in wine samples including 337 compounds already identified as wood markers (Supplemental Figure 8). The higher presence of wine-based compounds in whiskies older than 12 years, could correspond to finishing practices that could take place after 10 years of ageing independently of the region of production.

Conclusion

In conclusion, FT-ICR-MS based metabolomics was used as a powerful tool to discriminate whiskies according to maturation time and geographical origin, regardless of the manufacturer, the barrel used and the production process. Untargeted metabolomics allowed to study a broad range of metabolites and the extraction of markers that were significantly different depending on the geographical origin and/or the maturation time. A prediction model was further developed for the comprehensive evaluation of maturation time authenticity. Interestingly, the flow of transfer of compounds from wood barrels to distillates is not constant and homogeneous over the maturation time. The largest transfer of compounds from the barrel to the whisky was observed between ten and twelve years of maturation with up to 1182 molecular formulas that showed a significant increase in their peak intensities as a function of the maturation time. 467 of these were identified as wood-specific compounds and 338 are typically also found in wine. The explanatory hypothesis for the higher presence of wine-related compounds in whiskies older than 12 years, compared to whiskies younger than 12 years could be the finishing practice that takes place after 12 years of ageing independently of the region of production.

Figure 1 .

 1 Figure 1. Several thousands of signals were found between m/z 100 to 600. The enlargements at m/z 261 exemplarily show the differences in chemical composition between the four geographical areas (Figure 1). Between m/z 261.0 and 261.2, whisky spectra only differed in the intensities of the eight detected peaks. For example, signals at m/z 261.04044 and at m/z 261.17075 correspond to the molecularformulas [C13H9O6] -and [C13H25O5] -, respectively,

Figure 1 :

 1 Figure 1: Visualization of ESI(-) FT-ICR-MS spectra of four samples from Highlands, Lowlands, Speyside and Islay A) in the mass range m/z 100 -600. B) The enlargement at the nominal mass m/z 261 shows eight ions detected and assigned to unique molecular formulas. C) The Venn diagram shows the number of detected features which are common and specific to all samples analyzed for the four regions. D) Van Krevelen diagram (H/C vs. O/C) of the common composition found in all analyzed whisky samples. Points in van Krevelen diagrams are colored according to their elemental composition: CHO in blue, CHOS in green, CHON in orange, and CHONS in red. Scaling is according to the mean peak intensity.

Figure 2 :

 2 Figure 2: Principal component score plots of 106 whisky samples from Scotland, analyzed by FT-ICR-MS and colored according to A) the geographical origin and B) the maturation time (from 6 up to 18 years). The first two components explained 38.4% of the total variability.

Figure 3 :

 3 Figure 3: Map of the four whisky regions in Scotland (Highland, Speyside, Lowland and Islay). PLS regression model based on 106 whisky samples and van Krevelen diagrams showing the extracted "VIP" features for each region, respectively.

Figure 4 :

 4 Figure 4: A) Visualization of ESI (-) FT-ICR-MS spectra (m/z 100 to 600) of four samples matured in barrel for 6, 10, 14 and 18 years. B) Enlargement of the nominal mass at m/z 261 showing eight masses detected and assigned to unique molecular formulas. C) PLS model colored according to maturation time, based on 71 whisky samples (R²Y(cum)=0.57 and Q²(cum)=0.53). Comparison of prediction results versus the maturation time for D) the samples used to build the PLS model and for E) the 5 predicted samples. Van Krevelen diagram (H/C vs. O/C) of the significant features, which contribute to the discrimination of maturation time including F) features which showed an increase in peak intensity and G) features with decreasing peak intensity as a function of the maturation time, respectively. H) Bar chart of the absolute intensity of some biomarkers that showed a significant change in peak intensities. Points in the van Krevelen diagrams are colored according to their elemental composition, CHO in blue, CHOS in green, CHON in orange, and CHONS in red. Scaling is according to the mean peak intensity.

Figure 1 Figure 4

 14 Figure 1

Figure 2

Figure 3