Sébastien Giordano

Simon Bailly

Loic Landrieu

Nesrine Chehata

Temporal Structured Classification of Sentinel 1 and 2 Time Series for Crop Type Mapping

Keywords: classification, temporal regularization, conditional random fields, agriculture, monitoring, Sentinel images

As part of the EU Common Agricultural Policy (CAP) reform of 2020, each EU member country is expected to suggest new farmland management protocols. Currently, farmers must manually declare each year their crop types into the Land-Parcel Identification Systems (LPIS), a geographic information system identifying the land use of agricultural parcels within each EU member country. Such a protocol remains tedious and error-prone. Automatic Earth observation image analysis can help achieving such a task. Leveraging the recent availability of precise and frequent Sentinel acquisitions, this work aims to automate the LPIS update.

We propose modeling the crop type of parcels from a sequence of (radar and optical) satellite acquisitions, as well as LPIS entries of previous years, with a linear-chain Conditional Random Field. The novelty lies on the fusion of multi-modal images at the feature level and the integration of temporal knowledge extracted from existing land-cover databases. We tested our model on two large-scale French study areas (≥1250 km 2 ), which are geographically distant and show different agronomic rules: the Seine et Marne (North of France) and the Alpes de Haute-Provence (South East). We use a granular nomenclature comprised of 25 categories.

Our model demonstrates promising results for the task of automating the LPIS update: 89.0% overall accuracy is reached in Seine et Marne (10 categories of the 25 present on the area) and 72.9% in Alpes de Haute-Provence (14 categories). We show that the temporal modeling increases the accuracy by +2.6% and +4.6%, respectively.

Introduction

Motivation

The Sentinel 1&2 satellites provide open and free acquisitions exhibiting unprecedented characteristics which are perfectly tailored to agriculture monitoring. Most critically, the high temporal frequency (5-6 days) is very beneficial for identifying crop types. The Sentinel program will be maintained at until at least 2030, which allows us to chronicle both short and long-term evolutions. Besides, the multispectral Sentinel-2 images display relevant spectral bands to agricultural monitoring, that are complementarity to C-Band Sentinel-1 radar images. Finally, the high spatial resolution of both optical and radar images (10-20 m) authorizes parcel-level approaches. For all these reasons, the European Commission strongly recommended the use of the Sentinel programme for reshaping the procedures to monitor the Common Agricultural Policy (CAP) European Commission [START_REF]Towards future Copernicus service components in support to agriculture?[END_REF]. This therefore becomes a significant line of research in the forthcoming decade.

In Europe, several use-cases of agricultural monitoring using Sentinel images have been proposed [START_REF]Towards future Copernicus service components in support to agriculture?[END_REF], such as crop monitoring (crop area estimates, crop map products, crop phenology indicators), controlling CAP payments with remote sensing (permanent grasslands, greening measures, . . .), updating, and controlling the quality of the Land Parcel Identification System (LPIS) or precision farming at the farm-level.

In this paper, we focus on automating crop type mapping. In Europe, the agricultural land cover information is manually updated yearly by the farmers themselves. They input the type and surface area of their parcels. This manual declaration is complicated for farmers, error-prone and leads to expensive control procedures by external agencies. In this context, a pre-filled agricultural declaration based on supervised classification techniques would allow the farmers to only have to validate the declaration, cutting down on control and input errors. The automation of crop type declarations requires a robust classification model, based on Sentinel images observations but also on ancillary data such as LPIS archives to improve the classification results. Indeed, the crop type identification may be improved using a priori knowledge on management practices and especially on crop rotations per parcel. The choice of the crop type and agricultural practices on a parcel are strongly dependent on past events over previous years. The LPIS archives provide such information. Modeling these temporal structures in combination with the Sentinel image time series can lead to significant gains in classification accuracy. In this paper, we focus on benefiting from crop rotation knowledge.

Related work 1.2.1. Multi-temporal satellite images for crop mapping

In the literature, many studies showed the potential of multi-temporal satellite images for crop type mapping [START_REF] Palchowdhuri | Classification of multi-temporal spectral indices for crop type mapping: a case study in coalville, uk[END_REF]. Inglada et al. [START_REF] Inglada | Assessment of an operational system for crop type map production using high temporal and spatial resolution satellite optical imagery[END_REF] assessed stateof-the-art methods for automatic crop mapping with multi-temporal and very high spatial resolution optical images. For this purpose, five different classification approaches using SPOT4 and Landsat-8 images were compared on 12 different study areas worldwide. Best results (Overall accuracy of 80% for 6 annual crop classes) were obtained with the Random Forest classifier [START_REF] Breiman | Random forests[END_REF].

Immitzer et al. [START_REF] Immitzer | First experience with Sentinel-2 data for crop and tree species classifications in central Europe[END_REF] used mono-temporal Sentinel-2 images for agricultural and forest land cover classification. A multi-temporal approach has been proposed by Kussul et al. [START_REF] Kussul | Parcel-based crop classification in Ukraine using Landsat-8 data and Sentinel-1a data[END_REF]. Landsat-8 and Sentinel-1 time series were used on a study area in Ukraine. A pixel-based classification accompanied with a parcel-based regularization (majority voting) was proposed using LPIS ancillary data. An overall accuracy of 89% was reached but the nomenclature was limited to 6 annual crops (Winter wheat, Winter rapeseed, Maize, Sugar beet, Sunflower, Soybeans) and large parcels were generally considered (> 250 ha). More recently, experiments were led at the country level (Czech Republic) by the Sen2-Agri consortium [START_REF]Sen2-Agri, Czech agriculture national demonstrator -final report[END_REF]. A multi-sensor (Sentinel-1, Sentinel-2) pixel-based supervised classification was performed. The LPIS was used for both learning and validation steps. A crop map was produced every month using Sentinel-1 radar images (December to September) and Sentinel-2 optical images (March to September). The overall accuracy was greater than 80% and each land cover type had a F-score greater than 60%. The quality of the classification was further increased as more data was acquired. However, the nomenclature was here again very limited (7 classes) and does not fully integrate temporal knowledge from existing data.

Two major aspects of the state-of-art of crop type classification remain to be improved in order to obtain a reliable pre-filled declaration system: the precision on small parcels and the granularity of the nomenclature.

Crop rotation consideration

Crop rotations are known to be useful management practices, improving agricultural yields [START_REF] Berzsenyi | Effect of crop rotation and fertilisation on maize and wheat yields and yield stability in a long-term experiment[END_REF] and soil quality [START_REF] Karlen | Crop rotation effects on soil quality at three northern corn/soybean belt locations[END_REF]. References in literature on their role in agriculture are numerous [START_REF] Karlen | Crop rotations for the 21st century[END_REF]. To take into account crop rotations in crop type prediction, two questions have to be answered: (i) how to model the crop rotations? (ii) how to integrate crop rotations in a land cover classification process?

Modeling crop rotations

Two different approaches can be used to model rotations. The first one consists in using a priori agronomist expert knowledge. The second one would be to automatically learn crop rotations from the statistical analysis of past practices, as found in the LPIS archives. The crop rotation knowledge can then be modeled with a common representation in agronomy: transition matrices representing transition probabilities between crop types from year to year.

The ROTAT model [START_REF] Dogliotti | ROTAT, a tool for systematically generating crop rotations[END_REF] is based on expert knowledge. The tool automatically generates all possible crop rotations over an area and performs a selection taking expert knowledge into account. Castellazzi et al. [START_REF] Castellazzi | A systematic representation of crop rotations[END_REF] introduce a mathematical framework based on transition matrices to model crop rotations at the landscape scale. Dury et al. [START_REF] Dury | Models to support cropping plan and crop rotation decisions. a review[END_REF] review various crop rotation models and emphasize that these models are too static. To improve the models, the authors propose integrating different time scale dynamics. Indeed, crop rotation models have several limitations. The information is never available at the parcel scale. The models are strongly dependent on the study area; they cannot adapt to environmental or agricultural management changes. Some studies have been proposed to overcome these limitations. Aurbacher and Dabbert [START_REF] Aurbacher | Generating crop sequences in land-use models using maximum entropy and Markov chains[END_REF] underline the importance of having rotation information at the parcel-level and consequently propose to use Markov chains. Concerning adaptation capacities to different study areas, Detlefsen and Jensen [START_REF] Detlefsen | Modelling optimal crop sequences using network flows[END_REF] use network modeling techniques to find on a given parcel the optimal rotations knowing a selection of crops. Finally, the influence of environmental changes such as global warming on agricultural practices has been studied by Olesen et al. [START_REF] Olesen | Impacts and adaptation of European crop production systems to climate change[END_REF], then modeled by Aurbacher et al. [START_REF] Aurbacher | Influence of climate change on short term management of field crops-a modelling approach[END_REF]. Nevertheless, few studies take this parameter into account in crop rotation models.

Learning rotations using past data on crop types is a way of overcoming the limitations of a priori expert knowledge approaches.Both approches can also be successfully combined. For instance, the CarrotAge tool [START_REF] Ber | Studying crop sequences with CarrotAge, a HMM-based data mining software[END_REF] allows to perform a spatio-temporal analysis by training a Hidden Markov Models with a land cover database. Results are thereafter interpreted by an agronomist in order to be integrated in soil or water management studies. Another example is the CropRota method [START_REF] Schönhart | CropRota-a crop rotation model to support integrated land use assessments[END_REF] that combines agronomic criteria and land cover information to generate crop models at the farm or region levels. The ROTOR model [START_REF] Bachinger | ROTOR, a tool for generating and evaluating crop rotations for organic farming systems[END_REF] relies on control sampled farms, field surveys and expert knowledge.

These three models can take into account changes in agricultural practices or study areas. However, expert knowledge is still required, and such models can not be used at the parcel level.

The parcel-level crop information provides a finer crop type nomenclature but needs heavy field surveys and more complex models than those described above. Xiao et al. [START_REF] Xiao | Modeling the spatial distribution of crop sequences at a large regional scale using land-cover survey data: A case from France[END_REF] used such parcel-based information to describe the spatial distribution of crop sequences at a large regional scale, but with a limited 3-class nomenclature. For many years, several European countries as well as the United States maintain and annually update a geographical information system on agricultural parcels, i.e. the Land Parcel Identification System (LPIS) [START_REF] Boryan | Monitoring us agriculture: the us department of agriculture, national agricultural statistics service, cropland data layer program[END_REF]. Parcel-based crop rotation information at a national scale can be derived by exploiting past annual versions of the LPIS. The information on environmental or agricultural changes is then implicitly contained in these data. Several studies proposed to use the LPIS to study and model crop rotations. Leteinturier et al. [START_REF] Leteinturier | Adaptation of a crop sequence indicator based on a land parcel management system[END_REF] used several years of LPIS to compute crop rotation indicators. More recently, Osman et al. [START_REF] Osman | Assessment of a Markov logic model of crop rotations for early crop mapping[END_REF] produced and assessed a land cover classification based on the former version of the French LPIS.

Integrating crop rotations into classification pipelines

Only few studies have focused on the integration of crop rotation information into classification pipelines. Osman et al. [START_REF] Osman | Assessment of a Markov logic model of crop rotations for early crop mapping[END_REF] studied early crop mapping using only the LPIS. A prediction model based on Markov logic was proposed, but not in combination with remote sensing observations. Results showed that this model has better classification precision than models based on remote sensing observations at the beginning of the crop growing season. This is mainly due to the fact that few satellite images are available at the beginning of growing season and that crops are hard to distinguish at this development stage: observation-based classification is unreliable compared to temporal prediction. This statement remains to be confirmed with the use of higher temporal resolution Sentinel images. Other studies proposed to introduce a temporal structure, using the Hidden Markov Chains, in a classification pipeline but they aim at modeling phenology instead of crop rotations [START_REF] Aurdal | Use of hidden Markov models and phenology for multitemporal satellite image classification: Applications to mountain vegetation classification[END_REF][START_REF] Leite | Hidden Markov Models for crop recognition in remote sensing image sequences[END_REF][START_REF] Siachalou | A hidden Markov models approach for crop classification: Linking crop phenology to time series of multi-sensor remote sensing data[END_REF]. Modeling phenology is useful to detect optimal temporal intervals to acquire remote sensing images for crop mapping [START_REF] Conrad | Derivation of temporal windows for accurate crop discrimination in heterogeneous croplands of Uzbekistan using multitemporal RapidEye images[END_REF]. Kenduiywoa et al. [START_REF] Kenduiywoa | Spatial-temporal Conditional Random Fields crop classification from Terrasar-X images[END_REF] modeled phenology information into a Conditional Random Field, but the classification is performed at different dates through the year.

Specific objectives

The general objective of this article is to assess the feasibility of a pre-filled agricultural declaration based on Sentinel 1+2 images for crop types and surfaces. This paper focuses on crop type predictions. Related work described in Section 1.2, allow us to formulate the following specific objectives:

• Thematic contribution. Sentinel images are part of a recent satellite program (October 2014 and December 2015 for radar and optical images, respectively). Few research studies have yet assessed the crop type prediction accuracies that could be obtained with these time series. This study will especially focus on using an exhaustive national crop nomenclature, the consideration of all agricultural parcels (with no minimal crop area consideration), and the complementarity between Sentinel-1 and Sentinel-2 images.

• Methodological contribution. Crop rotations have rarely been combined with satellite observations into a classification pipeline. Our objective is to propose a method that integrates the crop rotation temporal structure into the classification process and to assess the capacity of the proposed method to improve crop type prediction accuracies.

Sentinel data, LPIS, and the study sites will be described in Section 2 as well as the necessary image and vector preprocessing steps. Two study sites with very different agricultural management practices and parcel sizes are chosen. In Section 3, we propose a parcel-based classification with a Random Forest classifier [START_REF] Breiman | Random forests[END_REF] and a temporal-structured framework to integrate crop rotation information. Results are given and discussed in Section 4. The models learned for the LPIS crop rotations and land cover prediction from the Sentinel images are first assessed independently. A combination of both models is finally evaluated and discussed.

Site and material

Study sites

In order to assess the feasibility of a pre-filled declaration of crop types, two complementary large-scale sites were chosen in France. The characteristics and location of each site are provided in Figure 1 and Table 1. site is characterized by a flat relief, with a large parcellar and a majority of cereal crops. This site is a permanent observatory of a Group of Scientific Interest; GIS Oracle (http://gisoracle.irstea.fr/). This site will be noted Site77 in reference to the national number of the corresponding region.

Table 1 illustrates the difference of both sites in terms of covered area and cultivated crop types. Figure 3 gives for both sites the normalized histograms of parcels area. Figures 4 and5, show Site04 and Site77, respectively, with the corresponding parcels and nomenclatures.

Land Parcellar Identification system

In France, the Land Parcellar Identification System is called Registre Parcellaire Graphique (RPG). It is available on the whole territory since 2002.

It corresponds to a Geographic Information System (GIS) of cultivated and non cultivated areas (NCA) that may correspond to isolated trees, hedges, groves, artificial areas . . . For cultivated areas, the RPG gathers the geometric information (i.e., the parcel delimitation) and the corresponding semantic information allowing to identify each agricultural parcel such as the owner, the operator, the area, the crop type etc (cf. Figure 2). The geometrically stable blocks of parcels from 2010 to 2014 (available LPIS data) were used to learn crop rotations. The latter step was processed on the corresponding departments of both sites ((77): 5915 km 2 and (04): 6925 km 2 ), in order to get more robust crop type transitions.

Nomenclature

The 28-class nomenclature used in the RPG declarations is driven agricultural monitoring needs: hence, it exhibits a fine-grained categorization.

However, some classes are indistinguishable from purely remote sensing obser-

Multimodal Sentinel-1 & Sentinel-2 images

Our framework is fed with multimodal and temporal Sentinel images dedicated to agricultural applications and environmental monitoring. The high spectral resolution of the images and high temporal sampling rate make these acquisitions particularly well-suited for crop mapping (Table 2).

Sentinel-1 (S1) sensor provides band-C SAR images while Sentinel-2 (S2) provides multispectral images.

Figure 6 shows a Sentinel-2 optical and a Sentinel-1 radar acquisitions on Site04, more precisely on Oraison Village.

Sentinel-2 (S2) is a multispectral sensor with 13 bands covering the VIS-SWIR domain, which measures the reflectance of surface objects in different optical domains. In particular, its near infra-red (NIR) and red-edge bands allow a fine characterization of crops. Sentinel-1 is a C-band SAR. This allows to measure scattering coefficients that are related to an emitted waveform (λ = 5.5 cm in C-band). The recorded energy depends on the characteristics of the encountered object (slope, roughness, humidity,etc.) and on the emitted waveform (wavelength, polarization). Sentinel-1 has different acquisition modes that differ according to observed surface, sensor spatial resolution and polarization. The available mode on the studied sites was the Interferometric Wide (IW) mode. This mode presents a dual polarization:

• VV : Incident waveform is polarized vertically and the antenna records it vertically. This polarization allows us to characterize the soil roughness.

• VH : Incident waveform is polarized vertically and the antenna records it horizontally. This polarization provides volumetric information on vegetation.

In addition, we use directly the GRD (Ground Range Detected) image format, which corresponds to the average of approximately five Single Look Complex ((SLC) acquisitions corrected by the incidence angle and resampled at 10 m spatial resolution (Section 2.5.2).

Sentinel image pre-processing

Figure 7 illustrates the optical and radar pre-processing steps to obtain parcel-based features that will be fed into our classification workflow.

Sentinel images repositories

Sentinel images are available on several platforms: the Copernicus scihub (https://scihub.copernicus.eu/) at the European level, the scihub mirror Peps (https://peps.cnes.fr) and the downstream service Theia (https: //www.theia-land.fr/) at the French national level. Sentinel-2 images (10 bands) were automatically downloaded from the Theia platform as they were is located where ascendant Sentinel-1 images overlap making the available radar images more numerous.

Sentinel 1 pre-process

The dual polarization GRD S1 images were first calibrated to σ 0 radar backscattering coefficient. Then, the orthorectification was performed using the SRTM digital terrain model and the georeferencing information supplied with the GRD files. The speckle is partly removed using a simple 5 × 5

Lee filter [START_REF] Lee | Digital image enhancement and noise filtering by use of local statistics[END_REF]. The information was then averaged at the parcel level. In addition to VV and VH radar features, an extra radar feature (

σ 0 V H σ 0 V V ) was
derived. This ratio is known to be more robust to acquisition system errors or environmental factors such as soil moisture. As a result, Veloso et al. [START_REF] Veloso | Understanding the temporal behavior of crops using Sentinel-1 and Sentinel-2-like data for agricultural applications[END_REF] argue that it is a more temporal stable indicator than the σ 0

V V et σ 0 V H backscattering coefficient.
Average and standard deviation of these three features per radar image are then computed for each date and for each parcel of the study sites. The number of obtained features are shown in Table 3 .

Site Nb of dates Optical features Radar features Total

Sentinel 2 pre-process

Sentinel-2 images downloaded from Theia platform were already orthorectified in a cartographic system and calibrated in TOC reflectance. When the study area covered more than one tile, the corresponding tiles and cloud masks were merged. On Site77, only 12 Sentinel-2 optical images were obtained in 2016 as shown in Figure 8 with corresponding cloud cover whereas 23 images are available on Site04. The missing data (clouds) were filled using a multi-temporal spline interpolation [START_REF] Inglada | OTB Gapfilling, a temporal gapfilling for image time series library[END_REF]. Average and standard deviation of the 10 spectral bands per optical image were then computed for each date and for each parcel of the study sites. The number of obtained features are shown in Table 3.

J F M A M J J A S O N D 74%1% 10% 99% 22% 0% 0% 43% 99% 4%98% 98%

LPIS pre-process

Seven editions of the French LPIS were used. From 2010 to 2014, the crop type is defined at the block level (Figure 2) (i.e., the majority crop type of the block) with a nomenclature of 28 classes. For 2015 and 2016 editions, the crop type is known at the parcel-level with a detailed nomenclature (over 300 classes). In order to use all available editions, all LPIS were aggregated to the 25 considered classes. Geometrically stable parcels were identified using a GIS spatial join between 2015 and 2016 versions. The number of stable parcels for both sites is given in Table 1. The 2016 stable parcels were used for training and validation of the supervised classification model (Section 3.1).

The 2015 stable parcels were necessary to train the temporal structured method (Section 3.2). Finally, a similar GIS spatial join was performed at the block level for 2010-2014 LPIS editions, over the whole corresponding regions (04 and 77) in order to get more robust crop transitions. Thus, a 5-year series of crop type transitions were obtained on stable blocks over both 04 and 77 regions. The number of transitions for each study site is given in Table 1. These transitions will be used in the temporal structure estimation (Section 3.2.1).

Methodology

Our method can be broken down into two components: parcel-wise classification and temporal modeling. The first part aims to predict the crop types per parcel using Sentinel time series with no temporal structure. The second component consists in integrating into a structured model the temporal structure derived from previous LPIS editions and based on crop rotations.

Parcel-based multi-source classification

The first step in our pipeline is crafting discriminative features from sequential satellite acquisitions. For each parcel, we consider all available (and interpolated) optical and radar data over one year. To obtain parcelbased rather than pixel-based features, we consider the average and standard deviation of each spectral band over the pixels composing the parcel. Those attributes are aggregated into a tensor of dimensions equal to twice the number of acquisitions over the year for each parcel.

We then perform feature selection by iteratively removing the features with the least importance until the cross-validated classification score starts decreasing over our learning set. For a given parcel i and a given year t we denote X (t)

i ∈ R D the tensor of aggregated selected features, with D the selected feature size. As certain classes were over-represented in our data sets, each class is weighted proportionally to the square root of the inverse of its number of instances.

A Random Forest classifier is selected for the classification task. It provides parcel-wise prediction under the form of a pseudo-probability.

Temporal-structured classification of parcels

We now consider the temporal structure of each parcel independently. We denote by

X i ∈ R T ×D the sequence of features X (t) i ∈ R d for the parcel i for the years t = 1, • • • , T .
Likewise, we denote Y ∈ K N ×T the labels of each parcel for each observed year with K the set of all possible labels.

Temporal structure

The aim of this step is to model the yearly crop rotations in order to improve crop type prediction.

As stated in 1.2.2, crop rotation has a significant impact on land cover. This dependency is modeled with a discriminative linear chain Conditional Random Field (CRF) of order m, as shown in Figure 9. For a parcel i, we model the posterior distribution P (Y i | X i ) of the labels Y i given the observed features X i as:

P (Y i | X i ) = 1 Z exp T t=1 O(Y (t) i , X i ) + T t=m+1 I(Y (t-m) i , • • • , Y (t) i , X) , ( 1
)
where Z is a normalizing factor, O the observation potential, and I the interaction potential, described below.

Observation potential: The observation potential models the link between the observed features and the label of each parcel. O(Y (t) i , X i ) is taken as the logarithm of the pseudo-probability given by the random forest classifier, described in Section 3.1. Interaction potential: This potential models the temporal dependencies between the parcel's labels. We model this potential as the transition probability from a sequence

Y (t-m) i , • • • , Y (t-1) i to a label Y (t)
i . For the sake of simplicity, we choose an homogeneous parameterization, independent of the observed features, and shared by all parcels and years:

I(Y (t-m) i , • • • , Y (t) i , X) = log M (Y (t-m) i , • • • , Y (t) i )) , with M ∈ R k m + a tensor such that {i 1 ,••• ,i m-1 }∈K m-1 M i 1 ,••• ,i m-1 ,im = 1 for all i m ∈ K.
This tensor can be interpreted as a transition probability from a sequence in K m-1 to a label in K [START_REF] Liu | Using local transition probability models in Markov random fields for forest change detection[END_REF].

Y (1) i ... Y (t-2) i Y (t-1) i Y (t) i X (1) i ... X (t-2) i X (t-1) i X (t) i

Learning

The observation potential is obtained by training the random forest classifier. Learning the transition tensor M from labeled data can be done in close form. For i 1 , . . . , i m ∈ K m , let:

Mi 1 ,...,im = N i 1 ,...,im N i 1 ,...,i m-1 ,
with N i 1 ,...,im the number of sequences i 1 , . . . , i m observed in the labeled data for all parcels and all years, and N i 1 ,...,i m-1 the number of sequences i 1 , . . . , i m-1

observed for the T -m first years.

To account for the large size of this matrix (k m ) and to prevent numeric issues, we perform a Laplacian smoothing with α = 1 as described in Manning et al. [34, 11.3.2].

Inference

The aim of this step is to predict the label Y (t) i from the labels and observations of the previous years for an unseen parcel. This can be directly computed by injecting the observation and interaction potentials obtained from the trained models into p(Y

(t) i = k | Y (t-m,••• ,t-1) i , X i ) ∝ exp O(k, X (t) i ) + I(Y (t-m) i , • • • , Y (t-1) i , k) ,
and normalizing the results to obtain a probability.

Results and discussions

In this section, we present the experimental setup and the evaluation metrics. Results are illustrated on both sites and will be discussed with regard to methodological and thematic objectives presented in Section 1.3.

Experimental setup

The random forest classifier is composed of 100 decision trees. The meta-parameters of the forest, such as the maximum number of attributes considered at each node, are chosen by k-fold cross-validation with k = 4.

For the temporal structure, spatio-temporal homogeneity hypothesis allows us to estimate the transition tensor M . For each study site, only the geometrically-stable parcel-blocks over the corresponding departments and a period of 5 years are used. The number of 5-year transitions that contribute to estimate M is given in Table 1.

The data is randomly split equally into a training and a testing sets. In practice, the availability of numerous LPIS editions allows us to have a high number of training parcels. The model is trained and validated on the training set while the quality of the model is estimated on the testing set. The overall accuracy (OA) is used to assess the general performance of the model, while the F-score for each class allows us to estimate the per-class quality. To sum up this information, the F-scores are averaged, with and without weighting by the class cardinality. The results are given as an average over 10 runs of the random forest classifier.

Transition matrix assessment

The transitions were computed on geometrically stable parcels between 2010-2014. In this study, only first order transitions were computed, i.e., between one given year and the previous year. Figures 10 and11 give the estimated transitions between crop classes as Hinton diagrams, for both sites.

On Site04, the most probable transitions are to and from permanent crops, such as olive groves, vignards, orchards, estive landes, meadows and fruit trees reaching 98.34%, 93.87%, 92.72%, 98.31%, 91.89% and 84.23%, respectively. From the transition matrix, we can observe that the standard rotation patterns of annual crops are generally not applied in this area. The rape seed, proteins and sun flowers have probabilities of 76.53%, 66.78% and 64.25%, respectively to be transformed to other cereals the following year.

On Site77, the most probable transitions are conversely observed for the annual crops. Agronomical rules for annual crops rotation seems to be much more enforced in this area. The rape seed and proteins have probabilities well-known 3-year rotation for farmers of this area. Permanent crops such as meadows, orchards, fruit trees, vineyards have a probability of being carried over the next year of 94.45%, 83.65%, 81.39% and 62.12%, respectively.

Quantitative and qualitative analysis

This section illustrates prediction accuracies using different feature combinations and both unstructured and structured approaches.

Quantitative evaluation on Site04

From the results displayed in The results confirm that optical and radar combination leads to the best results in both unstructured and structured approaches.However, the global results remain low for this site, with F-scores varying between 0.58 and 0.64. This is due to the very small size of parcels and highly imbalanced classes.

For example, parcels of permanent classes (meadows, fruit trees, vignards, olive groves) are overrepresented, as shown in Figure 4.

The confusion matrix for unstructured optical/radar configuration (Table 6) shows that most ambiguities occur on meadows classes, other cereals and forage crops. When integrating the temporal structure, the overall accuracy is improved by 4.6%. Moreover, while the OA using radar data is very low (63.9%), it is improved by the structured approach by 11.8%, confirming the impact of temporal structure even if the accuracy of the parcel-wise prediction is low.

In Table 5, we display the F-score, recall, and accuracy measures per class for both unstructured and structured approaches using radar and optical combination.

As for per-class accuracies, temporal structure highly improves the accuracies of permanent crops (fruit trees +36%, vignards +29.7%, olive groves +25.6%, aromatic groves +19%) reaching F-scores higher than 92.9%. The F-score for meadows (permanent and temporary) is improved by 4.1%. These However, F-scores of annual crops classes decrease when using temporal structure; only slightly so for corn, other cereal, rape seed and sunflowers but barley, other oilseeds, protein and forage crops are significantly more misclassified (cf. table 6). This is due to the fact that the transitions of annual crops are less stable and highly vary with regard to agricultural practices and operators in this area.

Impact of parcel size on Site04

Figure 12 shows the impact of parcel size on parcel-wise global accuracies.

Since Site04 is highly fragmented, when keeping only large parcels (area >3 ha), the number of parcels is reduced by 77.5%. Overall accuracies are highly improved by 15%, 5.9% and 8.7% for radar, optical and aggregated optical/radar attributes, respectively). Indeed, due to the limited spatial resolution of Sentinel-1 images, radar attributes are more sensitive to small radar attributes still give the best results for the unstructured approach. As with Site04, Table 8 shows that the meadows class is often confused with other classes, particularly fruit trees and other cereals. In this case, this is due to the combination of trees and bare soil found in meadows, which has a low volumetric radar response.

The structured approach improves overall accuracies by 2.7%, 4.6% and 2.6% for radar, optical, and aggregated optical/radar attributes, respectively.

The improvement of temporal structure is lower than for Site04, as the initial parcel-wise accuracies are already strong (OA>82%). However, one can observe that the structured approach leads to the best results with radar attributes only.

As for Site04, the best improvements occur on permanent crops such as meadows and fruit trees (cf. Table 9). Moreover, the temporal structure improves the prediction of some annual crops such as other cereals, Rape seeds and proteins since they have a high first order transition probability to other cereals.

Impact of parcel size on Site77

Figure 13 shows the impact of parcel size on the global accuracies of the parcel-wise prediction. Since Site77 is less fragmented than Site04, keeping only large parcels (area > 3 ha), reduces the number of parcels by % 52,3%.

As with Site04, when considering larger parcels, the overall accuracy is greatly improved with radar attributes (by 7.3% reaching 97.1%), which confirms the high sensitivity of radar images to parcel size.

For optical attributes, the relation between parcel size and accuracy is less clear. This may be due to the data imbalance, with some classes that are more present in small parcels (0.5-1.5 ha) and which are well identified such as Rape seed or Protein. Thus, removing these small parcels may decrease the overall accuracy. The classes that are responsible for the higher OA observed between 0 and 0.75 ha are Fiber plants and Meadows where meadows parcels are much smaller than for other classes and that fiber plants parcels are homogeneous in terms of area (cf. Figure 5).

As for radar data sensitivity to limited parcel sizes, some improvements of our framework could be undertaken on radar pre-processing. Indeed, the GRD Sentinel-1 images were used as input images when each GRD radar pixel is already an average of 5 pixels of the SLC (Single Look Complex). Using SLC orthorectification accuracy of Sentinel-1 images and allows us to correct radar scattering coefficient from the effect of terrain slope (γ 0 calibration). As a result, the intra-class variability of signal could be reduced.

Conclusion and perspectives

This study demonstrated the efficiency of an automatic prediction of crop types using Sentinel 1 and 2 images and LPIS archives as a first step for an automatic filled pre-declaration intended for farmers.

From a thematic point of view, this paper demonstrates the efficiency of multi-temporal and multi-source Sentinel (optical and radar) images for crop type classification on two different sites. Very satisfactory discrimination results are obtained even with a large number of classes (>10). The remaining limitations are essentially the parcel size with regard to the spatial resolution of images and the complexity of the nomenclature. Results highlight the high sensitivity of radar data to small parcel sizes, especially on Site04. This issue can be reduced by a refining the preprocessing framework of radar data.

Moreover, the complexity of certain agronomic classes such as forage crops should be discussed with agronomists in order to design a proper hierarchical classification.

From a methodological point of view, we integrated the temporal structure by automatically modeling the crop rotations using prevision editions of LPIS.

This model appears to be very efficient, improving the global classification prediction of crop types. However, the impact of integrating temporal structure varies highly among classes. Although a positive impact is demonstrated on permanent crops using first order crop transitions, this impact is fairly limited or even detrimental for annual crops. Other transition orders should be investigated to confirm the interest of temporal structure for annual crops.

Finally, thanks to the large volume of available LPIS data that can be used as ground truth, and the free availability of numerous Sentinel images, deep learning methods should be tested to learn parcels' embeddings.

  The first one is located in South Eastern France, in the Alpes de Haute-Provence region, in the Durance river Valley. It is representative of a Mediterranean cultivated area. This site will be called Site04 in reference to the national number of the corresponding region. It covers 1050 km 2 and is characterized by a highly variable topography, a very fragmented parcellar while giving a high diversity of crop types. The site is an observatory of the French mapping agency (IGN) where crop observations are made annually. The second one is located near Paris, in the Seine et Marne region, at the north of Coulommiers town. It covers 233 km 2 . Contrarily to the Site04, this

Figure 1 :

 1 Figure 1: Localisation of Site04 and Site77

Figure 2 :

 2 Figure 2: Details of the French LPIS (RPG). (a) Non cultivated areas (NCA); (b) a block of contiguous parcels, belonging to the same operator (dashed lines); (c) a Parcel, defined by its owner, operator, surface, crop type. . .

Figure 3 :

 3 Figure 3: Normalized histogram of parcel area for Site04 and Site77.

Figure 4 :Figure 5 :3

 45 Figure 4: Site04 : 2016 RPG parcellar superimposed to a very high resolution Digital Terrain Model.

Figure 6 :

 6 Figure 6: Site04, focus on Oraison village: Sentinel-1 (5 Aug. 2015). Composed color image with different polarizations related to retrodiffusion coefficients: σ vv , σ vh , σ vh σvv . Site04-Oraison village: Color-infrared Sentinel-2 (3 Aug. 2015). See text for more details.

Figure 7 :

 7 Figure 7: Sentinel-1 and 2 pre-processing steps

Figure 8 :

 8 Figure 8: S2 optical images over the year 2016 and corresponding cloud cover on Site77.

Figure 9 :

 9 Figure 9: Graph structure of the temporal dependency at order 2.

Figure 10 :

 10 Figure 10: Site 77 -Representation of the transition matrix with a Hinton diagram.

of 97 .Figure 11 :

 9711 Figure 11: Site 04 -Representation of the transition matrix with a Hinton diagram.

  images directly would allow us to achieve greater spatial resolution. Then,

Figure 13 :

 13 Figure 13: Unstructured parcel-based classification for year 2016, Site77 -Impact of parcel size on accuracy; x-axis: only parcels whose surface area exceeds the threshold (in ha) are considered.

Table 1 :

 1 Comparison of both study sites in terms of area and crop types.

	Class	# parcels -Site04 # parcels -Site77
	Corn	147	350
	Barley	517	158
	Other cereals	2176	889
	Rape seed	154	85
	Sun flower	293
	Other oilseeds	116
	Protein(peas)	87	76
	Fiber plants		76
	Forage crops	1215	46
	Meadows	3652	725
	Fruit trees	298	30
	Vignards	249
	Olive groves	1029
	Aromatic crops	1452
	Vegetables	520	131
	Total nb classes	14	10
	Total nb parcels	11905	2566
	Site area (km 2 )	1050	233
	Total nb stable parcels 2015/2016	9230	1902
	Total transitions 2010-14	29478	36891

vations, such as lands left in fallow for different amount of times. Consequently, we present in Table

1

a slightly adjusted nomenclature of 25 classes. 14 and 10 of these classes are present on Site04 and Site77, respectively.

Figure

3

shows the distribution of parcel sizes on both sites. One can see that Site04 is much more fragmented with very small parcels while Site 77 shows large parcels reaching 20 ha.

Figures

4 and 5

show the 2016 RPG edition i.e., the ground truth data on Site04 and Site77 respectively, with the corresponding classes on each site. One can observe that, on Site04, dominant crops are: cereals (23.8%), meadows (30.7%), aromatic crops (12.2%), forage crops (10.2%) and olive groves (8.6%). For Site77, two dominant crops are present: cereals (57.7%), meadows (28.3%), followed by vegetables (5.1%). The data are much more imbalanced in the latter case, making the 77 classification task more complex.

Table 2 :

 2 Characteristics of Sentinel sensors and used images: Sentinel-1, Sentinel-2.

Table 3 :

 3 Characteristics of the parcel-based features for both sites. See text for more details.

Table 4

 4

	, one can see that for parcel-based

Table 4 :

 4 Global prediction accuracies on Site04, using different feature and methodological configurations.

		Unstructured
	Config	OA F-score Weighted F-score
	Radar	0.639 0.587	0.610
	Optical	0.730 0.669	0.714
	Radar Optical 0.729 0.679	0.713
		Structured
	Radar	0.757 0.592	0.698
	Optical	0.776 0.627	0.725
	Radar Optical 0.775 0.645	0.723

Table 5 :

 5 Effect of temporal structure on Classification metrics on Site04, using aggregated radar and optical attributes.

		Unstructured		Structured
	Class	F score Recall Precision F score Recall Precision
	Corn	0.888	0.953	0.832	0.780	0.836	0.732
	Barley	0.393	0.848	0.256	0.185	0.673	0.107
	Other cereals	0.846	0.830	0.862	0.783	0.730	0.845
	Rape seed	0.923	1	0.857	0.712	1	0.554
	Sun flower	0.789	0.785	0.793	0.713	0.749	0.681
	Other oilseeds	0.570	0.745	0.463	0.408	1	0.258
	Protein	0.484	1	0.321	0	0.0714	0
	Forage crops	0.469	0.735	0.344	0.096	0.884	0.051
	Meadows	0.759	0.675	0.867	0.800	0.674	0.983
	Fruit trees	0.609	0.863	0.471	0.969	0.962	0.975
	Vignards	0.690	0.740	0.646	0.987	0.996	0.979
	Olive groves	0.737	0.819	0.670	0.993	0.997	0.989
	Aromatic crops 0.739	0.646	0.864	0.929	0.938	0.921
	Vegetables	0.611	0.725	0.529	0.674	0.891	0.542
	results were expected since the permanent crops have the highest transition
	probability as shown in Section 4.2.

Table 7 :

 7 Global prediction accuracies on Site77, using different attribute and methodological configurations.

		Unstructured
	Config	OA F-score Weighted F-score
	Radar	0.892 0.734	0.878
	Optical	0.824 0.624	0.809
	Radar Optical 0.890 0.744	0.885
		Structured
	Radar	0.919 0.776	0.911
	Optical	0.870 0.675	0.853
	Radar Optical 0.916 0.762	0.906

When looking at per-class accuracies in Table

8

, one can see that unlike Site04, very good overall accuracies are obtained for annual crops (Corn (94.1%), Barley (89.8%), Other cereals (94.7%), Rape seed (95.9%), . . . ) with the exception of forage crops. Indeed this class is hard to classify using satellite imagery since it is an agronomic class, making it harder to identify using only spectral or radar scattering information even on temporal images.

Table 8 :

 8 Effect of temporal structure on Classification metrics on Site77, using aggregated radar and optical attributes.

	Class	Unstructured		Structured
	Class	F-score Recall Precision F-score Recall Precision
	Corn	0.941	0.929	0.953	0.878	0.831	0.935
	Barley	0.898	0.937	0.862	0.816	0.785	0.849
	Other cereals 0.947	0.956	0.9378	0.954	0.941	0.968
	Rape seed	0.959	0.975	0.944	0.969	0.985	0.954
	Protein	0.949	0.932	0.968	0.953	0.967	0.939
	Fiber plants	0.974	1	0.950	0	0.1	0
	Forage crops	0	0.1	0	0.705	0.778	0.648
	Meadows	0.868	0.814	0.930	0.955	0.943	0.967
	Fruit trees	0.010	0.090	0	0.940	1	0.892
	Vegetables	0.895	0.914	0.877	0.451	0.970	0.297

Acknowledgments

The authors want to thank the ASP (French payment agency) and the French agricultural ministry for the rich discussions and their guidance for the nomenclature choice and crop rotation schemes.

Quantitative evaluation on Site77

Overall accuracies and F-scores are displayed in Table 7. They are averaged over 10 runs of the random forest classifier. Contrarily to the previous site, on Site77, radar attributes give better results than optical attributes (OA = 0.892 Vs. 0.824) for the unstructured approach. This is due to a combination of frequent acquisition problems and high cloud cover in 2016, leading to many missing optical Sentinel-2 data on this study area (cf. Figure 8). In addition, the parcels on Site77 are larger and thus more compatible with radar Sentinel-1 image spatial resolution. However, aggregated optical and

simple speckle filter was used [START_REF] Lee | Digital image enhancement and noise filtering by use of local statistics[END_REF] on a restricted local neighborhood (5×5). This is not a problem for large parcels as the radar scattering coefficients are averaged afterwards at the agricultural parcel level. However, when parcel surface areas approach the Sentinel-1 spatial resolution, this method is no longer suitable. The robustness of adaptive radar speckle filtering to small objects should be investigated [START_REF] Deledalle | NL-SAR: A unified nonlocal framework for resolution-preserving (pol)(in) SAR denoising[END_REF]. Finally, using a very high spatial resolution Digital Terrain Model instead of the SRTM could improve absolute


Muriel Gasc-Barbier

Véronique Merrien-Soukatchoff

GASC-BARBIER Muriel -description of Ax-les-Termes solpe cutting for a Rock fall hazard benchmark case Description of Ax-les-Thermes cuttings for a Rock fall hazard Benchmark case

Keywords: benchmark, rock slope, natural hazard

This paper will present the site, the geological and structural data that constitute a benchmark on rock fall hazard. All the data provided come from a site that was studied for the bypass of Ax-les-Thermes (Ariège, France).

AIM OF THE BENCHMARK

Starting from fracturing data and a slope (or underground) geometry, estimates of the rock fall hazard can be ran under different methodologies, using various software. The idea of this benchmark is to compare practices in:

 Regrouping fractures in families:  Consideration or not of sampling bias  Mechanical properties taken into account for the discontinuities (deterministic or stochastic)  Rock fall hazard estimation  The stability indicators given from the stability analyses. These indicators can be varied for example the results can be expressed in terms of:

-biggest unstable block -more frequent block (possibly with a confidence interval) -total unstable volume (weight) and dispersion around this volume -number of unstable blocks -blocks displacement and movement types -local or Global Safety Factor  The inclusion of a reinforcement (bolting) and the presence of water in the hazard assessment  software used for these different steps  How the geometrical, and mechanical properties uncertainties are taken into account

We propose to work on data collected before excavating a high cutting in a migmatitic rock excavated for a bypass in the French Pyrenees where 861 fracture measurements were performed [START_REF] Gasc-Barbier | Design of large cuttings in jointed rock[END_REF]. Another benchmark on rock fall hazard in underground context will also be proposed based on the data collected before and during the excavation of St Beat tunnel, also in the French Pyrenees.

DATA SURVEY

Location

Ax-les-Thermes lies 130 km south of Toulouse, just 30 km from the border between France and Andorra, in the primary axial zone of the French Pyrenees. Three orogenic phases that affected the entire mountain lead in considerable fracturing of the site, which is apparent in different forms depending on the area studied.

Geological information

The bypass encounters Quaternary superficial formations, fluvio-glacial materials and then the rock substratum, where the benchmark takes place. The slope under study (Esquiroulet cuttings) is mainly composed of augen gneiss and leptinic gneiss, in practice it is impossible to separate them into two distinct areas and hence the formation is best described as migmatite (in the widest sense). All these rocks are heavily fractured. GASC-BARBIER Muriel -description of Ax-les-Termes solpe cutting for a Rock fall hazard benchmark case

Scanlines and borehole images

For several reasons, excavation of the site was stopped when the proposed 42 m high cutting was 15 m high. This provided an opportunity to realize scanlines on the new cuttings to obtain statistical information on the joint properties. A scanline is a horizontal line (of about 10 m) drawn on a rock face [START_REF] Priest | Discontinuity analysis for rock engineering[END_REF]. The line can, in fact, be vertical or horizontal but should have spatial references (dip and dip direction if possible). As recommended by many authors [START_REF] Aftès | Recommandations relatives à la caractérisation des massifs rocheux utile à l'étude et à la réalisation des ouvrages souterrains[END_REF][START_REF] Priest | Discontinuity analysis for rock engineering[END_REF][START_REF] Khanlari | Instability assessment of slopes in heavily jointed limestone rocks[END_REF], the standard properties of all the joints that intersect the line are mapped: dip, dip direction, spacing, trace length, weathering, aperture, infill and roughness. Several lines in different directions were drawn in order to try to intersect the maximum number of joint sets. In this particular case and because of the geometry of the site, the scanlines were drawn on three faces as depicted in Figure 1 but their lengths were quite different. It was not possible to use vertical lines because of the difficulty of access. In total, 11 scanlines (from 6 to 30 m) were plotted; including 856 discontinuities over 181.5 m. Borehole images, also allow having a kind of "scanlines". Two boreholes were inspected by camera allowing an enhancing of the data and "internal" data of the rock mass. All the data can be found on https://www.researchgate.net/profile/Muriel_Gasc-Barbier/publications.

Other provided information

Other data are provided on https://www.researchgate.net/profile/Muriel_Gasc-Barbier/publications: Elastic waves velocities, uniaxial and Brazilian tests and shear tests on natural joints.

Other provided information

Figure 1a shows a 3D view of the site when scanlines were drawn. Figure 1b gives the proposed geometry of the cutting to be reinforced for the benchmark.

Figure 1 :

 1 Figure 1: Ax-les-Thermes cuttings -a-localization of the different scanlines, b-proposition of geometry for the benchmark

GASC-BARBIER Muriel, Cerema, Toulouse, F-31400, muriel.gasc@cerema.fr

MERRIEN-SOUKATCHOFF Véronique, Cnam, Paris, F-75002, veronique.merrien@lecnam.net


Alice Gadea

Anne-Cécile Le Lamer

Sophie Le Gall

Catherine Jonard

Solenn

Ferron

Daniel Catheline

Damien Ertz

Pierre Le Pogam

Joël Boustie

Françoise Lohezic -Le Devehat

Maryvonne Charrier
email: maryvonne.charrier@univ-rennes1.fr

Françoise Lohezic

INTRATHALLINE METABOLITE PROFILES IN THE LICHEN Argopsis friesiana SHAPE 1

Keywords: Herbivory, chemical ecology, lichen-gastropod interactions, Subantarctic islands, 38 Stereocaulaceae, Notodiscus hookeri. 39 40

Lichen-gastropod interactions generally focus on the potential deterrent or toxic role of 24 secondary metabolites. To better understand lichen-gastropod interactions, a controlled feeding 25 experiment was designed to identify the parts of the lichen Argopsis friesiana consumed by the 26 Subantarctic land snail Notodiscus hookeri. Besides profiling secondary metabolites in various 27 lichen parts (apothecia, cephalodia, phyllocladia and fungal axis of the pseudopodetium), we 28 investigated potentially beneficial resources that snails can utilize from the lichen 29 (carbohydrates, amino acids, fatty acids, polysaccharides and total nitrogen). Notodiscus hookeri 30 preferred cephalodia and algal layers, which had high contents of carbohydrates, nitrogen, or 31 both. Apothecia were avoided, perhaps due to their low contents of sugars and polyols. Although 32 pseudopodetia were characterized by high content of arabitol, they were also rich in medullary 33 secondary compounds, which may explain why they were not consumed. Thus, the balance 34 between nutrients (particularly nitrogen and polyols) and secondary metabolites appears to play 35 a key role in the feeding preferences of this snail.

INTRODUCTION 43

Lichens are above all a lifestyle generally presented as a stable organism corresponding to the 44 symbiotic association of fungi with photobionts [START_REF] Gargas | Multiple origins of lichen symbioses in fungi 495 suggested by SSU rDNA phylogeny[END_REF]. Fungal symbioses with algal 45 or cyanobacterial photobionts are called chloro-and cyanolichens, respectively. Sometimes the 46 mycobiont coexists with both green algae and cyanobacteria and forms tripartite lichens, 47 referred to as cephalolichens, because cyanobacteria are located inside specialized fungal 48 compartments (cephalodia). 49

In a lichen thallus, green algae produce polyols while cyanobacteria synthesize sugars. These 50 photosynthetically produced carbohydrates may be metabolized by the fungus partner into 51 mannitol and/or arabitol (Richardson et al., 1968). The mycobiont also produces secondary 52 metabolites (or specialized metabolites) that accumulate in the cortex or in the medullar layer 53 [START_REF] Molnár | Current results on biological activities of lichen secondary 535 metabolites: a review[END_REF]. In most chlorolichens, the acetyl-polymalonyl pathway leads to the 54 synthesis of secondary metabolites, mainly phenolic compounds such as depsides and 55 depsidones [START_REF] Stocker-Wörgötter | Metabolic diversity of lichen-forming ascomycetous fungi: culturing, 584 polyketide and shikimatemetabolite production, and PKS genes[END_REF]. In most cyanolichens and in cephalodia containing 56 cyanobacteria, terpenoids and mycosporines predominate [START_REF] Rundel | The ecological role of secondary lichen substances[END_REF][START_REF] Roullier | Characterization and 560 identification of mycosporines-like compounds in cyanolichens. Isolation of mycosporine 561 hydroxyglutamicol from Nephroma laevigatum Ach[END_REF]. 57

The co-occurrence of two photobionts in cephalolichens is reflected by the diversity and 58 allocation of metabolites produced, cephalodia being devoid of phenolic compounds, which 59 are found in the other parts of the thallus, but rich in N-containing compounds (Roullier et al. 60 2011). As a general trend, a higher content of nitrogen occurs in cephalo-and cyanolichens as 61 the atmospheric nitrogen fixed by the cyanobacteria is converted to amino acids by the 62 mycobiont [START_REF] Rai | Cyanolichens[END_REF]. 63

Phenolic compounds in lichens have various roles. Some may deter lichen feeders such as 64 molluscs [START_REF] Coker | Damage to lichens by gastropods[END_REF]Asplund et al. 2010b), mammals [START_REF] Cook | Paresis and death in elk (Cervus elaphus) due 481 to lichen intoxication in Wyoming[END_REF][START_REF] Nybakken | Lichen compounds restrain lichen 544 feeding by bank voles (Myodes glareolus)[END_REF]) 65 and arthropods [START_REF] Gerson | Lichen-arthropod associations[END_REF][START_REF] Wieners | Carbon-based secondary compounds in the lichen 593 Hypogymnia physodes deter detrivorous woodlice[END_REF]. Various physical and chemical factors may 66 influence the attractiveness of lichens for lichenivorous organisms, such as roughness and 67 growth form, occurrence of defence compounds and N-fixing ability of the lichen (Asplund and 68 Wardle 2016). In studies focussing on interactions between plants and phytophagous animals, 69 different hypotheses have been proposed postulating a trade-off between the cost of operating 70 chemical defence and the benefit of enhancing survival. In his work on plant-herbivore 71 coevolution, [START_REF] Mckey | Adaptive patterns in alkaloid physiology[END_REF] developed the "optimal defence theory" (ODT) which was later 72 applied by Asplund et al. (2010b) to lichen-gastropod interactions. This theory predicts that 73 secondary metabolites in plants are allocated proportionally to the risk of attack for a specific 74 tissue and the value of this tissue in terms of the reduction in fitness that its loss would entail. In 75 many lichens, secondary metabolites occur in higher concentrations in reproductive parts, 76 whether sexual (apothecia) or asexual (isidia, soralia), than in the rest of the thallus, or occur 77 only in these reproductive parts. [START_REF] Hyvärinen | Phenols in reproductive and somatic 516 structures of lichens: a case of optimal defence?[END_REF]. [START_REF] Rundel | The ecological role of secondary lichen substances[END_REF] conceptualized another 78 hypothesis, suggesting that the higher the nitrogen content, the higher is the lichen's palatability. 79

Hence, lichens with low nitrogen content and high concentration of secondary metabolites 80 (particularly depsides and depsidones) should suffer low grazing. Rundel's hypothesis was 81 confirmed for most chlorolichens, which were indeed more grazed after non-destructive removal 82 of the secondary metabolites by the acetone-rinsing method [START_REF] Solhaug | Acetone rinsing-a method for testing ecological and physiological 577 roles of secondary compounds in living lichens[END_REF]83 Asplund and Wardle 2013). Conversely, Rundel's theory was not validated for some 84 cyanolichens which were not consumed although characterised by high nitrogen levels and low 85 content of secondary chemicals. The presence of intracellular repulsive compounds or of 86 metabolites not extractable by acetone (i.e. mycosporines) was hypothesized, due to the 87 persistence of a deterrent effect after acetone rinsing [START_REF] Solhaug | Secondary lichen compounds as protection against excess solar 579 radiation and herbivores[END_REF]. 88

Lichens are widespread from cold to tropical habitats, from wet to xeric conditions and in 89 nutrient-poor habitats, such as alpine and polar areas, where they are often dominant (Sanders 90 2001). In Subantarctic islands, where lichens are abundant, it could be important for 91 phytophages to exploit such an abundant resource and thus allow an expansion into new 92 environments. The Subantarctic land snail Notodiscus hookeri, an exclusive lichen feeder, is 93 particularly widespread on Crozet (Possession Island) and Kerguelen Archipelagos (Charrier et 94 al. 2013). This species provides an appropriate model organism for studying lichen-gastropod 95 interactions. Recently, [START_REF] Gadea | Which specialized metabolites does the native 489 subantarctic gastropod Notodiscus hookeri extract from the consumption of the lichens 490 Usnea taylorii and Pseudocyphellaria crocata?[END_REF] demonstrated the generalist status of N. hookeri by a 96 comparative study of snail metabolism after feeding on Usnea taylorii (chlorolichen) and 97

Pseudocyphellaria crocata (cyanolichen). These authors observed several snail strategies to 98 overcome potentially toxic metabolites present in lichens. 99

Optimal foraging in generalist species, such as most terrestrial gastropods, requires a balance 100 between the access to useful nutrients and the need to overcome toxic compounds. Therefore, 101 the snail should be able to select lichen parts according to their nutritional quality and to avoid 102 toxic compounds. In this study, we focussed on the cephalolichen Argopsis friesiana Müll. Arg. 103 (Stereocaulaceae), a tufty ≤ 4 cm tall species, widespread in low-nutrient, windy fell-fields of 104 Subantarctic islands [START_REF] Lamb | The lichen genus Argopsis Th[END_REF][START_REF] Galloway | The lichen genera Argopsis and Stereocaulon in New Zealand[END_REF]. Argopsis friesiana is characterized by a rigid 105 pseudopodetium, cylindric-corraloid phyllocladia, concave-scutelliform apothecia and botryose 106 cephalodia. Its main photobiont is a green alga (Trebouxia s. lat.), its secondary autotrophic 107 partner is the cyanobacterium Stigonema sp., located in the cephalodia. To understand how this 108 snail could benefit from eating this lichen, we examined the influence of lichen chemistry on snail 109 nutrition. For this purpose, we explored how the metabolites were distributed in the different 110 parts of the lichen A. friesiana and how this distribution shaped the snail grazing patterns. 111

Intrathalline sectorization of the symbiotic partners facilitated the separation of A. friesiana into 112 four different parts to be analysed separately. Three hypotheses were tested: The snail (1) 113 avoids lichen parts with high concentrations of secondary metabolites, (2) prefers lichen parts 114 that are rich in essential nutrients (sugars, amino acids, fatty acids, polysaccharides) and (3) 115 prefers thallus parts rich in nitrogen. 116 117

MATERIALS AND METHODS 118

Indian Ocean, was collected on Possession Island (Crozet archipelago) at « Le Donjon » 120 (46°25'01"S, 51°50'15"E, 200 m, DON) during the austral summer 2015. The identification of the 121 species was confirmed by one of us (DE) by comparison with type material and the species was 122 sequenced for the first time. A sequence of nuclear ribosomal internal transcribed spacer (ITS) 123 of the material used in this study was deposited in GenBank (MG947382) as a reference for our 124 material because ITS is used as a universal DNA barcode marker for fungi [START_REF] Schoch | Nuclear ribosomal internal transcribed spacer 571 (ITS) region as a universal DNA barcode marker for Fungi[END_REF]. 125

The species resembles the sympatric Stereocaulon cymosum Cromb., which differs by its 126 transversely septate ascospores versus muriform ascospores in A. friesiana. Voucher 127 specimens were deposited at the herbarium of the Faculty of Pharmacy of Rennes 1, 128

Department of Pharmacognosy and Mycology, under the reference REN000140. 129

To investigate the palatability of A. friesiana, snails were collected from sites for which 130 we had an official collection permit. Seventy-five adults (shell size ≥ 5 mm) were collected during 131 2015 austral summer at Mascarin (46°26'10.09"S; 51°45'20.58"E, 600 m -930 m) (Online 132

Resource 1, Fig. S1). The specimens, found in fell-fields including the lichen habitats, belong to 133 the soft-shell ecophenotype described by [START_REF] Charrier | Soil calcium availability influences shell 467 ecophenotype formation in the Sub-Antarctic land snail, Notodiscus hookeri[END_REF]. 134 135 Lichen Preparation. Argopsis friesiana was cut into four parts as explained in Fig. 1a. The 136 pseudopodetium consists of (i) a central axis, the medulla, made up of compact fungus hyphae 137 (Pf) and (ii) an outer layer with fungi and green algae (Pa). Phyllocladia are the terminal 138 branchlets of pseudopodetia containing algae (Ph). Cyanobacteria live inside numerous black 139 cephalodia (Ce), external gall-like structures emerging laterally from the pseudopodetium. 140

Apothecia, the sexual reproductive parts of the lichen, are scattered on the upper surface of the 141 thallus. The underside of the apothecium is similar to the algal layer of the pseudopodetium. For 142 chemical analyses, algal layer and phyllocladia were combined because of their structural 143 similarity. 144 145 Experimental Design. Snails were reared in boxes (7.5 x 6 x 4.5 cm) under a thermal cycle, 6°C 146 during night (10h) and 10°C during day (14h). A moist synthetic foam in the bottom (5 mm thick) 147 ensured snail activity. After four days of starvation, the snails (n=75) were divided into 15 groups 148 of five individuals. They were allowed to feed upon A. friesiana for two nights (48 h experiment). 149

Each group received a fertile thallus, on which apothecia and cephalodia were counted. Snail 150 consumption was controlled after one night and damage was measured after two nights under a 151 stereomicroscope (Stemi 2000-C, Zeiss, France). 152

Level of grazing (LD) was adapted from [START_REF] Nimis | Lichen chemistry and selective grazing by the coleopteran Lasioderma 542 serricorne[END_REF] and measured as 153 follows: LD0 intact, LD1= low grazing (<25%), LD2=medium damage (>25-50% consumed) and 154 LD3= high damage (>50%). For pseudopodetia, grazing on the algae-containing cortical layer 155 (Pa) was assessed separately from that on the fungus axis (Pf) and the levels of grazing were 156 defined as: LD0 = algal layer and fungal axis of the pseudopodetia remained intact, LD1 = 157 minimum damage (algal layer consumed <25%, fungal axis intact), LD2 = medium damage 158 (algal layer consumed <50%, fungal axis intact), LD3 = high damage (algal layer Pa and fungal 159 axis Pf consumed, > 50%). To assess the damage to cephalodia and apothecia, a grazing ratio 160 was determined as follows: Number of apothecia or cephalodia with grazing marks / Number of 161 apothecia or cephalodia available. To compare grazing on cephalodia and apothecia with 162 damage to algal layer, fungal axis and phyllocladia, the levels described above (LD in %) were 163 analyses of each part of the lichen were made using the ambient air source 196 enabling the chemical profiling of unprocessed pieces of thallus (Le [START_REF] Pogam | Analysis of the cyanolichen Lichina 527 pygmaea metabolites using in situ DART-MS: from detection to thermochemistry of 528 mycosporine serinol[END_REF][START_REF] Pogam | In situ DART-MS as a versatile and rapid 524 dereplication tool in lichenology: chemical fingerprinting of Ophioparma ventosa[END_REF]; 197 more detailed protocol in Online Resource 2). The identification of metabolites was done against 198 compounds formerly identified within the genus Argopsis and the family Stereocaulaceae (Bodo 199 and Molho 1974;[START_REF] Lamb | The lichen genus Argopsis Th[END_REF][START_REF] Galloway | The lichen genera Argopsis and Stereocaulon in New Zealand[END_REF]. 200 201 Primary Metabolite Profiling. Soluble carbohydrates, polyols and amino acids were extracted 202 from each part of the lichen, following the methodology of [START_REF] Gravot | Diurnal oscillations of metabolite abundances 502 and gene analysis provide new insights into central metabolic processes of the brown 503 alga Ectocarpus siliculosus[END_REF]. Twenty mM of 203 adonitol (soluble carbohydrates and polyols) and 10mM of 3-aminobutyric acid (BABA) (amino 204 acids) were used as internal standards. Soluble carbohydrates and polyols were profiled by a 205

Gas Chromatography -Flame Ionization Detector (GC-FID) (Thermo-Fisher Scientific, Waltham, 206 CA, USA) according to [START_REF] Adams | Simultaneous determination by capillary 434 gas chromatography of organic acids, sugars, and sugar alcohols in plant tissue extracts 435 as their trimethylsilyl derivatives[END_REF] and [START_REF] Lugan | Metabolome and water status phenotyping of 530 Arabidopsis under abiotic stress cues reveals new insight into ESK1 function[END_REF]. Amino acids were detected 207 with an Ultra High Performance Liquid Chromatography -Diode Array Detector (UPLC-DAD, 208 Waters), by using an AccQ-Tag Ultra derivatization kit (Waters Corporation, Milford, MA, USA) 209 following the protocol described in [START_REF] Gravot | Diurnal oscillations of metabolite abundances 502 and gene analysis provide new insights into central metabolic processes of the brown 503 alga Ectocarpus siliculosus[END_REF]. The metabolites were identified from 210 their retention time by comparison with external standards. Concentrations (mg.g -1 DM) were 211 obtained on the basis of internal standards. 212 213 Fatty Acids Profiling. Samples were extracted three times in a parallel synthesis apparatus 214 (Heidolph Parallel Synthesis Apparatus®-Instruction Synthesis 1, France), for 1 h each, with 215 chloroform/methanol (2:1, v/v, 6 mL) at 80°C, under constant shaking at 700 rpm (Vu et al. 216 2016). To quantify fatty acids (FAs), 100 µg of heneicosanoic acid (C21:0) were added as an 217 internal standard. The organic phase containing the total lipids was recovered and evaporated to 218 dryness under a nitrogen flow. Sample derivatization and GC-FID/GC-MS (Agilent 6890N, 219 Toulouse, France) analyses were conducted using the protocol of [START_REF] Vu | Gas chromatographic analysis to compare the fatty 590 acid composition of fifteen lichen species, with a focus on Stereocaulon species[END_REF]. Identification 220 of fatty acid methyl esters (FAMEs) was based on the retention times. The quantity of total FAs 221 (in mg.g -1 DM) was calculated using the internal standard. Identification and quantification of polysaccharide neutral sugars were performed by GC 232 after sulfuric acid degradation [START_REF] Hoebler | Rapid acid hydrolysis of plant cell wall 509 polysaccharides and simplified quantitative determination of their neutral 510 monosaccharides by gas-liquid chromatography[END_REF]. Five mg of AIR were dispersed in 13 M 233 sulfuric acid for 30 min at 30°C, then hydrolysed in 1 M sulfuric acid (2 h, 100°C). Sugars were 234 converted to alditol acetates according to [START_REF] Blakeney | A simple and rapid preparation of alditol 462 acetates for monosaccharide analysis[END_REF] and chromatographed on a TG-235 225 GC Column (30 x 0.32 mm ID) using TRACE™ Ultra Gas Chromatograph (Thermo 236 Scientific™; Waltham, CA, USA; temperature 205°C, carrier gas H2). Inositol was used as 237 internal standard. The quantity of each sugar (in mg.g -1 DM) was determined, using the 238 calibration curves drawn from standard solutions. Apothecia were completely avoided by N. hookeri (Fig. 1a-2). Similarly, the fungal axis 264 was not consumed, while the cortical algal layer and phyllocladia were grazed (Fig. 1a-2). The 265 levels of damage (LD) on cortical algal layer and phyllocladia ranged from < 25% to < 50% and 266 were similar (Fisher's exact test, P=0.710) (Fig. 2). Cephalodia were the preferred part, more 267 than 40% of them being highly damaged (= LD3) and significant grazing differences being 268 observed between cyanobacterial (cephalodia) and green algal parts (cortical algal layer and 269 phyllocladia) of the lichen (Fisher's exact tests, P ≤ 0.001) (Fig. 1b-2, Fig. 2). The mean (± s.e.) 270 number of cephalodia grazed reached 4.9 ± 0.5 in 48 hours and was dependent on the total 271 number of cephalodia per thallus (Spearman's rank correlation, rho = 0.59, P = 0.021). Three main compounds were identified: the depside atranorin (m/z 373) and two chlorinated 276 compounds (Table 1, Fig. S2). These two chlorodepsidones were identified as the dichlorinated 277 argopsin (m/z 395) and as the trichlorinated caloploicin (m/z 401) through isolation and 278 spectroscopic data analysis (Purification protocol and NMR data are given in Online Resource 2, 279 Fig. S3-S4, Table S1-S2). In situ DART-HRMS analysis confirmed identification of chlorinated 280 compounds detected by LC-DAD-MS from unprocessed pieces of lichens through exact mass 281 measurements (Online Resource 2, Table S3). Algal layers had the highest concentration of 282 secondary metabolites, particularly for atranorin, reaching 10.6 ± 2.3 mg.g -1 DM. In comparison, 283 atranorin varied from 0.5 ± 0.2 mg.g -1 DM in cephalodia to 1.1 ± 0.6 mg.g -1 DM in apothecia. In 284 pseudopodetia, the medullar caloploicin with argopsin accounted for about two-thirds of the 285 secondary metabolites, while atranorin alone approximately reached this proportion in apothecia 286 and in green algae-containing tissues (Table S4). 287 288 Among all detectable metabolites, 36 were found in apothecia, 33 in the medullar axis of 289 pseudopodetia, 41 in phyllocladia and algal layer of pseudopodetia and 42 in cephalodia (Table 290 1). The multivariate analysis (PPLS-DA) showed partitioning in the metabolite profiles among 291 parts of the lichen A. friesiana (Fig. 3a, Fig. 3b, contributions to PC1 axis, Online Resource 3 292 Fig. S5). The cross model validation used for the discrimination of the lichen parts (8 293 components, 999 permutations, NMC = 0.042; P = 0.001) indicated statistically different 294 chemical profiles. The PC2 axis distinguished the presence or the absence of photobionts while 295 the biont identity was discriminated on the first axis (PC1). Green algal layers were 296 characterized by proline and GABA, and by the secondary metabolites atranorin, argopsin and 297 caloploicin. On the opposite side of the first axis, discrimination of the cephalodia was mainly 298 due to trehalose, FAunsat (n-7) family, galactose, mannose and rhamnose from polysaccharides, 299 amino acids such as lysine and arginine, and by high nitrogen concentration (TN). The position 300 of cephalodia was negatively correlated with the abundance in glucose issued from the 301 degradation of polysaccharides. On the second axis (contribution to PC2 axis, Online Resource 302 Fig. S6), the position of apothecia was negatively associated with the occurrence of amino acids 303 and polyols (especially arabitol) and a positive but limited association with some FA (FAsat and 304 FAunsat (n-9) families) was observed. The position of the fungal axis of pseudopodetia on the third 305 axis (PC3) was positively correlated with the presence of fatty acids (FAsat, C18:2 n-6, C16:1 n-306 9), mannitol) and α-alanine (Fig. 3c, Fig. 3d, contributions to PC3 axis, Online Resource 3 Fig. 307 S7). 308 309 DISCUSSION 310 311 Among the three major secondary metabolites detected in Argopsis friesiana, the depside 312 atranorin and the chlorodepsidone argopsin have already been reported in this species (Huneck 313 and Lamb 1975). [START_REF] Yosioka | The structure of caloploicin, a new lichen trichloro-596 depsidone[END_REF] revealed the chlorodepsidone caloploicin in the genus 314

Caloplaca, but here, for the first time, caloploicin has been isolated from the genus Argopsis with 315 the complete NMR assignment for this compound (Online resource 2). 316

When consuming A. friesiana, N. hookeri at first meets external layers high in secondary 317 metabolites. Despite the abundance of atranorin in the cortex, it had no deterrent effect on snail 318 grazing. In phytophagous arthropods, atranorin was reported to elicit a significant avoidance 319 [START_REF] Nimis | Lichen chemistry and selective grazing by the coleopteran Lasioderma 542 serricorne[END_REF] or to reduce growth of insect larvae [START_REF] Slansky | Effect of the lichen chemicals atranorin and vulpinic acid upon feeding and 574 growth of larvae of the yellow-striped armyworm, Spodoptera ornithogalli[END_REF][START_REF] Pöykkö | Removal of lichen secondary metabolites affects food 547 choice and survival of lichenivorous moth larvae[END_REF]) 320 but no protective effect of atranorin was noticed in snails [START_REF] Hesbacher | Sequestration of lichen compounds by three 507 species of terrestrial snails[END_REF][START_REF] Gauslaa | Lichen palatability depends on investments in herbivore defence[END_REF]. 321

In the soft bodies of various snail species, [START_REF] Hesbacher | Sequestration of lichen compounds by three 507 species of terrestrial snails[END_REF] detected a compound 322 presumably originating from hydrolysis of atranorin and suggested its sequestration in snail 323 tissues. While the retention of lichen compounds in bodies of snails was not analysed as part of 324 the current study, this hypothesis is consistent with the previously reported ability of N. hookeri 325 to hydrolyse ester bonds of the tridepside tenuiorin, due to the alkalinity of its gut lumen (Gadea 326 et al. 2017). According to [START_REF] Lawrey | Chemical defense in lichen symbioses[END_REF], cortical compounds (i.e. atranorin in A. friesiana) are 327 assigned a protective role against high irradiance, whereas a number of medullary secondary 328 metabolites deter lichenivores. As caloploicin and argopsin are depsidones, they were expected 329 to occur predominantly in the medulla. In A. friesiana, much higher concentrations of these 330 compounds were found in the external green algal layers, but expressed in percentages they 331 represented less than 25% of the detected secondary metabolites in this part. By contrast, 332 caloploicin and argopsin reached 70% of the secondary metabolites found in the central axis of 333 the pseudopodetia. The medulla of A. friesiana remained intact after snail feeding, consistent 334 with a deterrent role of these compounds. This is in agreement with [START_REF] Asplund | Snails avoid the medulla of Lobaria pulmonaria and L. scrobiculata due to 440 presence of secondary compounds[END_REF], who 335 demonstrated that snails avoid the medulla of Lobaria species due to the presence of 336

depsidones. 337

According to the Optimal Defence Theory [START_REF] Mckey | Adaptive patterns in alkaloid physiology[END_REF]Asplund et al. 2010b), 338 reproductive tissues should be well defended. However, our results did not support the theory. In 339 A. friesiana, apothecia and cephalodia had similar concentrations of medullary secondary 340 metabolites. Unlike cephalodia, apothecia were avoided by N. hookeri, suggesting that (i) other 341 undetected toxic compounds deter snail grazing, or (ii) their avoidance may be due to low 342 nutrient content, or (iii) the roughness of apothecia may present a physical barrier to radula 343 motion. Therefore, the food choices made by N. hookeri cannot be explained simply by the 344 avoidance of secondary metabolites, supporting the conclusion of [START_REF] Lawrey | Vulpinic and pinastric acids as lichen antiherbivore compounds: Contrary 522 evidence[END_REF]. Other factors 345 or compounds such as essential nutrients for snail growth and reproduction might be involved. 346

Cephalodia and algal layers were selected by N. hookeri when given a choice among the 347 four parts of A. friesiana. These lichen parts were characterized by their richness in some 348 nutrients: (i) polyols and free sugars, among which trehalose was encountered in cephalodia 349 only, (ii) amino acids, among which proline was exclusive of algal layers. Avoidance of apothecia 350 by the snail might be linked to their low content of phagostimulant compounds. However, we 351 cannot explain why the snails did not consume the fungal axis of pseudopodetia. Although rich in 352 arabitol, the fungal axis was characterized by a lack of diversity in nutrients compared to that 353 found in apothecia. Moreover, the central axis made of compact hyphae, as well as apothecia, 354 may be resistant to grazing. Indeed, the radular teeth of N. hookeri are smooth and more or less 355 bulging or flat, features of non-specialized radulae potentially sensitive to rough (apothecia) and 356 fibrous (fungal axis) textures (Fig. S8). By contrast, the soft and gelatinous structure of the 357 cephalodia could promote grazing. This observation is in line with the study of [START_REF] Baur | Ultrastructure of snail grazing damage to calcicolous 460 lichens[END_REF], 358 who compared the ultrastructure of snail radulae and concluded that hook-shaped teeth, such as 359 those of Chondrina clienta, are necessary to scrape efficiently on crustose lichens. 360

We demonstrated the high palatability of cephalodia and algal layers for N. hookeri and 361 that the more numerous the cephalodia, the more consumed they were. Various factors could 362 explain these feeding choices. Due to their low levels of defences, cephalodia may be preferred 363 by snails [START_REF] Renner | Études microspectrophotométriques sur les céphalodies de Pseudocyphellaria 554 faveolata[END_REF]Asplund and Gauslaa 2010). Second, the snails could acquire a 364 balanced intake of essential nutrients by consuming several parts of the lichen with intrathalline 365 chemical specificities. If so, the snail would have to perceive nutrient quality even in the case of 366 sub-lethal toxicity particularly due to the occurrence of argopsin and caloploicin in A. friesiana. 367

An ability to select foods that contain essential organic or mineral nutrients despite ingestion of 368 secondary metabolites was evidenced in slugs by [START_REF] Cook | The influence of nutritional status on the 478 feeding behaviour of the field slug, Deroceras reticulatum (Müller)[END_REF] and in snails by Chevalier 369 et al. (2003). 370

The glucose residues in polysaccharides, although they were half as abundant in 371 cephalodia as in other parts, they suggested the presence of glucans, probably isolichenan (α-in lichen palatability following artificial N deposition, and [START_REF] Asplund | The impact of secondary compounds and functional 450 characteristics on lichen palatability and decomposition[END_REF] found that 399 N-enriched lichens were less palatable. However, these authors showed that the effect on 400 consumption rate depended on an interaction between nitrogen-fixing ability and growth form. 401

The cephalodiate and fruticose A. friesiana could be included in the group of fruticose N-fixing 402 species described by [START_REF] Asplund | The impact of secondary compounds and functional 450 characteristics on lichen palatability and decomposition[END_REF], which are consumed more than non-N-fixing 403 fruticose species. Hence, our results do not contradict those of [START_REF] Asplund | The impact of secondary compounds and functional 450 characteristics on lichen palatability and decomposition[END_REF]. 404

Because the snail tissues contain 10% nitrogen [START_REF] Speiser | Food and feeding behaviour[END_REF], we hypothesize that N. hookeri 405 should benefit from the higher N concentrations (1.6%) in cephalodia than in other lichen 406

tissues. 407

Studies on lichen-gastropod interactions have often examined only the repellent role of 408 secondary metabolites, but we showed here that the contents in nitrogen and primary 409 metabolites related compounds (such as amino acids, sugars and polyols) should not be 410 underestimated. Moreover, physical criteria such as the hardness of lichens thalli could be able 411 secondary metabolites, amino acids, free sugars and polyols, sugar-forming polysaccharides, 618 saturated and unsaturated fatty acids. Saturated fatty acids (FAsat) included C14:0, C15:0, 619 C16:0, isoC17:0, C18:0, C20:0, C22:0, C23:0, C24:0. Abbreviations used are given in Table 1. 620

For clarity, in the presence of collinearity (≥ 95%), only one of the two correlated variables was 621 considered (Table S3). c. Score plot 2-3 and d. corresponding loading plots. C14:0, C15:0, C16:0, isoC17:0, C18:0, C20:0, C22:0, C23:0, C24:0. Abbreviations used are given in Table 1. For clarity, in the presence of collinearity (≥ 95%), only one of the two correlated variables was considered (Table S3). c. Score plot 2-3 and d. corresponding loading plots.

  Chemical Analyses. To determine a sectorial distribution of metabolites, 166 three samples of several A. friesiana thalli (about one gram each) were prepared. Each thallus 167 was cut into the four parts apothecia, cephalodia, fungal axis and algal layers. The 12 samples 168 were ground separately with mortar and pestle in liquid nitrogen. To quantify secondary 169 metabolites, freeze-dried samples (40 mg each) were used. Primary metabolites, including 170 amino acids, free sugars and polyols, were analysed with 10 mg of each thallus part, while 50 171 mg of each part was needed for quantification of fatty acid and sugar-forming polysaccharides 172 and for determination of nitrogen content. 173 174 Secondary Metabolite Profiling. Each part of the air-dried lichen thalli was extracted with 100% 175 acetone (3 x 500 µL). Acetone extracts were solubilized in tetrahydrofuran (THF), filtered and 176 injected (10 µL at 0.5 mg.mL -1 ) in a Liquid Chromatography -Diode Array Detector (LC-DAD) 177 (Shimadzu®, Marne La Vallée, France, LC-20AD SP, injector SIL-20A HI, column oven CTO-178 20A, diode array detector SPD-M20A) and a mass spectrometer (MS) (Advion® expression 179 CMS, Ithaca, USA). Ionization was performed by negative-ion mode electrospray (ESI-). A 180 gradient system was applied: A (0.1% formic acid in water) and B (0.1% formic acid in 181 acetonitrile) for 48 min on a C18 column (Phenomenex®, Kinetex 2.6µ C18 100A, temperature 182 40°C) at a flow rate of 0.5 mL.min -1 : initial, 20% B; T 0-5 min, 20% B linear; T 30 min 80% B; T 183 35-42 min, 100% B linear; T 45-48 min, 20% B linear. Linear equations based on the calibration 184 curve were used to quantify atranorin and the results were expressed in mg.g -1 Dry Mass (DM). 185 Atranorin (JB/A/142) used for calibration belong to the library of single lichenic compounds of 186 UMR 6226, being isolated during previous phytochemical investigations. For the other 187 metabolites, also detected by UV on HPLC chromatograms, a relative quantification was 188 performed, by multiplying the area under the curve (AUC) and the extraction yield of the lichen 189 extract. Full scan mass spectra were recorded in negative-ion mode in a mass range of 100 to 190 1200 Da, applying the following parameters: detector gain 1200, ESI voltage 3.5 kV, capillary 191 voltage 180 V, source voltage 20 V, source voltage dynamic 20 V, nebulizer gas pressure 60 192 psig, desolvation flow gas rate 4 L.min -1 , capillary temperature 250°C and gas temperature 193 20°C. Data processing and evaluation for MS measurement were performed with the Data and 194 Mass Express 2.2.29.2 software (Advion). As a complement to LC-DAD-MS experiments, in situ 195

222 223

 222 Extraction of Sugar-Forming Polysaccharides. Approximately 50 mg of each dried lichen tissue 224 were reduced to a fine powder with a pestle and transferred in 22 mL cells of an accelerated 225 solvent extraction unit ASE® 350 (ASE™ 350, Thermo Scientific™ Dionex™), dried at 40°C 226 overnight under vacuum over P2O5, prior to weighing. Alcohol-insoluble residues (AIRs) were 227 obtained by using 80% ethanol at 2 mL.min -1 flow. The conditions for the ASE extraction were 228 set at 100°C for 20 min, followed by rinse volume of 150%, and a purge time (N2) of 30 s. Cells 229 of the ASE® 350 containing AIRs were dried at 40°C overnight under vacuum over P2O5, prior to 230 weighing, in order to determine the extraction yield. 231

  TN) Measurement. TN (in % of DM) was determined with an automated dry 241 combustion method (Dumas method ; Rhee 2001) using an elemental analyser (vario MICRO 242 cube, Elementar Analysensysteme, Hanau, Germany). 243 244 Statistical Analyses. A Fisher's exact test was applied to LD values obtained for damage to algal 245 layer and fungal axis of the pseudopodetia, phyllocladia, cephalodia and apothecia, in order to 246 compare the observed frequencies among lichen parts. In the text, mean ± standard deviation is 247 given. A Spearman correlation was used to compare LD values with the number of cephalodia or 248 apothecia carried out by the thallus. 249 To analyse the metabolic data, we first proceeded to fourth-root transformation of the 250 values, centring and autoscaling, because null values and great differences in concentrations 251 according to metabolite profiles were observed (van den Berg et al. 2006). Powered Partial 252 Least Squares -Discriminant Analysis (PPLS-DA) was performed on the normalized values to 253 discriminate the four parts of A. friesiana on the basis of their metabolic profiles (function cppls(), 254 package "pls" (Indahl et al. 2009)). The significance of the discrimination was tested using a 255 cross-model validation (functions MVA.cmv() and MVA.test()), developed in the package 256 "RVaidememoire" (Hervé 2016). 257 All statistical analyses were made using R software V. 3.4.0 (R Core Team 2017). 258 259 RESULTS 260 Snail Feeding Preferences. Snails from Mascarin site had shell sizes of 4.9 ± 0.5 mm (mean 261 ± s.d.). Shell size did not differ among groups (n = 15 x 5 individuals per group; ANOVA, F14,60 = 262 0.415, P = 0.964). 263

272 273Metabolite

 272 Profiling of A. friesiana. The chemical profiles in secondary metabolites of each part 274 of the lichen A. friesiana were established by LC-DAD-ESI-MS from their acetone extracts. 275

600 Fig. 1 Fig. 2 Fig. 3

 600123 FIGURES 599 600 Fig. 1 Argopsis friesiana (Stereocaulaceae) entire thallus (a) and cephalodia (b) before and after 601 snail consumption. Phyllocladia (Ph), algal layer of pseudopodetia (Pa) and cephalodia (Ce) 602 were consumed while apothecia (Ap) and fungal axis of pseudopodetia (Pf) remained 603 untouched. The specimens were dry before consumption (a-1 and b-1) and wet after snail 604

  data are related to this article: see Supplementary data file 626 Online resource 1. Fig. S1: Study sites on Possession Island, which belongs to Crozet 627 Archipelago (45° 30' -46° 30' S; 50° 00' -52° 30' E) in the Subantarctic region. 628 Online Resource 2: Chemical and spectroscopic data. Fig. S2: Three main secondary 629 metabolites of Argopsis friesiana: the chlorodepsidones argopsin (1) and caloploicin (2) and the 630

FiguresFigureFig. 1

 1 Figures INTRATHALLINE METABOLITE PROFILES IN THE LICHEN Argopsis friesiana SHAPEGASTROPOD GRAZING PATTERNS.

Fig. 2 Fig. 3

 23 Fig. 2 Differential grazing and level of damage (LD) to Argopsis friesiana by the snail Notodiscus hookeri. N = 15 lichen thalli given to snails. Pa = algae-containing layer of pseudopodetia, Ph = phyllocladia, Ce = cephalodia; apothecia and medullar part of pseudopodetia that remained untouched are not shown. The level damage classes are LD1 = minimum damage (>0% and <25%), LD2 = medium damage (25%-50%) and LD3 = high damage (>50%).

Table 1 Data set of the compounds identified in the four parts of the lichen Argopsis friesiana (Ap= apothecia, Ce= cephalodia, Ph= phyllocladia, Pa= algal layer of pseudopodetia and Pf= fungal axis of pseudopodetia): secondary metabolites, free sugars and polyols, amino acids, fatty acids and sugar forming polysaccharides. Values are means (n=3 replicates) in mg.g -1 dry mass with (minimum -maximum) in brackets, except the secondary metabolites, which are expressed in Area Under the Curve (samples: 0.5 mg.ml -1 ) *extraction yield (x10 5 ). Total polysaccharides (Total PS) and total nitrogen (TN) are expressed in percentage of DM. The sign '-' means that the metabolite was not detected. Galactose (PS_Gals) 63.48 (60.25 -68.95) 55.37 (48.63 -61.75) 56.35 (55.56 -57.29) 104.80 (89.72-116.90) Glucose (PS_Glc) 620.27 (580.91 -654.66) 718.34 (678.57 -785.53) 667.19 (630.97 -705.97) 316.64 (290.01 -337.77) Mannose (PS_Mans) 116.54 (112.02 -122.58) 77.24 (66.43 -85.07) 74.41 (73.73 -75.31) 182.82 (150.65 -201.89) Rhamnose (PS_Rhas) 4. 15 (3.89 -4.54) 4.02 (3.66 -4.48) 3.62 (2.88 -4.06) 6.53 (5.00 -7.29) Xylose (PS_Xyls) 9. 76 (7.93 -11.16) 7.00 (4.20 -10.39) 5.36 (4.99 -5.78) 10.04 (8.15 -12.60) Total PS (%) 72.39 (67.83 -76.95) 71.94 (66.17 -77.71) 66.55 (63.64 -69.46) 54.27 (50.60 -57.63) Total Nitrogen (TN) (%) 0.68 (0.57 -0.85) 0.37 (0.36 -0.38) 0.66 (0.63 -0.72) 1. 62 (1.46 -1.78)