Anh Khanh Cong

Yu Le

Frédéric Ogawa

Florent Dubreuil

Karim Grimaud

Gregory Mazeau

Shivalika Ziegler

Yoshiharu Tanwar

Gabrielle Nishiyama

Luc Potocki-Veronese

Choisnard

Le

Yu Ogawa

Frédéric Dubreuil

Florent Grimaud

Karim Mazeau

Gregory R Ziegler

Shivalika Tanwar

Yoshiharu Nishiyama

Gabrielle Potocki-Veronese

Luc Choisnard

Denis Wouessidjewe

Jean-Luc Putaux
email: jean-luc.putaux@cermav.cnrs.fr

Crystal and molecular structure of V-amylose complexed with ibuprofen

Keywords: amylose, ibuprofen, inclusion compound, crystal structure

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

Introduction

Amylose is a mostly linear homopolymer of α(1,4)-D-glucose that is extracted from native starch granules [START_REF] Buléon | Self-association and crystallization of amylose[END_REF]. Amylose-like linear polymers can also be synthesized in vitro by amylose synthases such as amylosucrase [START_REF] Potocki De Montalk | Amylosucrase from Neisseria polysaccharea: novel catalytic properties[END_REF] or starch phosphorylase [START_REF] Ohdan | Enzymatic synthesis of amylose[END_REF]. Amylose can be crystallized from solution by cooling or addition of a large variety of organic complexing agents [START_REF] Putseys | Amylose-inclusion complexes: Formation, identity and physico-chemical properties[END_REF][START_REF] Tomasik | Complexes of starch with organic guests[END_REF]. Depending on the molecular weight and concentration of amylose, nature and concentration of the complexing agent, and crystallization temperature, so-called "V-amylose" crystallosolvates with different morphologies and allomorphic structure are formed. In these crystals, the complexing molecules can be entrapped between amylose helices and/or inside their hydrophobic central cavity [START_REF] Tan | Starch-guest inclusion complexes: Formation, structure, and enzymatic digestion[END_REF]. This property can find applications, for instance, in the food industry (flavor encapsulation, texturing, etc.) or in pharmacology (vectorization, controlled release) (Carbinatto, Ribeiro, Colnago, Evangelista, & Cury, 2013;[START_REF] Conde-Petit | Structural features of starch-flavor complexation in food model systems[END_REF][START_REF] Kumar | Deciphering structures of inclusion complexes of amylose with natural phenolic amphiphiles[END_REF][START_REF] Ma | Formation of inclusion complexes of starch with fatty acid esters of bioactive compounds[END_REF].

β-cyclodextrins (β-CDs, cyclic oligo-amylose made of 7 α(1,4)-D-glucosyl units) have often been used to form inclusion compounds with various drugs in order to increase their solubility in water [START_REF] Eddine | Investigation of the inclusion complex of tolfenamic acid with β-cyclodextrin: Geometry and NBO analysis[END_REF]. In particular, complexes of β-CDs with ibuprofen (IBU), a widely used non-steroidal anti-inflammatory, analgesic and antipyretic agent, have been studied [START_REF] Mura | Interactions of ketoprofen and ibuprofen with β-cyclodextrins in solution and in the solid state[END_REF][START_REF] Salústio | Release profile of ibuprofen in β-cyclodextrin complexes from two different solid dosage forms[END_REF] and several methods to efficiently mix both constituents have been explored [START_REF] Khan | Ibuprofen β-cyclodextrin inclusion complexes: Evaluation of different complexation methods[END_REF][START_REF] Pereva | Efficiency of "cyclodextrinibuprofen" inclusion complex formation[END_REF]. It was shown that IBU molecules enter the hydrophobic cavity of the β-CDs, replacing the initial water molecules [START_REF] Pereva | Efficiency of "cyclodextrinibuprofen" inclusion complex formation[END_REF].

In line with the similitude with the hydrophobic inner cavity of cyclodextrins, the potential of the helical cavity of amylose to encapsulate poorly water-soluble bioactive molecules has also been considered [START_REF] Cai | Evaluation of amylose used as a drug delivery carrier[END_REF][START_REF] Feng | Evaluation of different analysis methods for the encapsulation efficiency of amylose inclusion compound[END_REF]Janaswamy, 2004). In particular, recent articles have reported on the preparation of inclusion compounds with IBU by acidification of an alkaline amylose solution (Yang et al., 2013;Zhang et al., 2016) or by in situ enzymatic polymerization (Yang et al., 2013). These studies agreed on the good in vitro stability of the complexes in a simulated gastric fluid and the gradual but sustained release of IBU by enzymatic degradation of the amylose matrix in a simulated intestinal medium (Yang et al., 2013;Zhang et al., 2016). Other studies have considered the entrapment of IBU in various types of starch matrices (Marinopolou, Papastergiadis, & Raphaelides, 2019;[START_REF] Mohammadi | Mesoporous starch aerogels production as drug delivery matrices: synthesis optimization, ibuprofen loading, and release property[END_REF]. However, while in some of these works, the X-ray diffraction profiles indicated the formation of a V-type semicrystalline complex, the crystallinity of the samples was very low. The substrates were mixtures of crystalline and amorphous regions in various amounts and it was not possible to estimate which fraction of the prepared substrates actually participated to the release experiments. The determination of a reliable molecular model supported by structural data collected from highly crystalline samples is mandatory in order to shed some light on the interaction of the guest molecules with the host amylose helices.

In this study, we have prepared model single crystals of VIBU complexes from dilute solutions of native amylose and fractions of a shorter amylose-like polymer enzymatically synthesized in vitro. The influence of parameters such as chain length, temperature, or amylose and ligand concentration was investigated. The morphology, structure and thermal properties of the lamellar single crystals have been characterized by complementary imaging, diffraction, spectroscopy and calorimetry techniques, and the data have been used to propose a plausible molecular model of the VIBU crystalline complex.

Experimental section

Materials

Native potato amylose was purchased from Sigma-Aldrich and purified as previously described [START_REF] Le | Polymorphism of crystalline complexes of amylose with fatty acids[END_REF]. Its weight-average degree of polymerization (𝐷𝑃 # $$$$$ ) was 2514 [START_REF] Le | Polymorphism of V-amylose co-crystallized with aliphatic diols[END_REF], so this fraction will be referred to as DP2500 in the following. Using the method summarized in the Supplementary Data section, a linear amylose-like polymer was enzymatically synthesized in vitro from sucrose by the amylosucrase from Neisseria polysaccharea, and subsequently fractionated [START_REF] Potocki-Veronese | Amylose synthesized in vitro by amylosucrase: Morphology, structure, and properties[END_REF]. We have selected fractions with 𝐷𝑃 # $$$$$ 60 and 130

(referred to as DP60 and DP130) (Table 1 and Supplementary Data Fig. S1). Amylose fractions with 𝐷𝑃 # $$$$$ 28, 192 and 601 (referred to as DP30, DP200 and DP600, respectively) synthesized in vitro by a glycoside-phosphorylase were purchased from Shoko Co. Ltd (Japan). The characteristics of all fractions are summarized in a synthesized in vitro by a glucoside-phosphorylase [START_REF] Ohdan | Enzymatic synthesis of amylose[END_REF] b synthesized in vitro by the amylosucrase from Neisseria polysaccharea [START_REF] Potocki-Veronese | Amylose synthesized in vitro by amylosucrase: Morphology, structure, and properties[END_REF] c native potato amylose characterized as described in [START_REF] Le | Polymorphism of V-amylose co-crystallized with aliphatic diols[END_REF] d values provided by the manufacturer

Preparation of crystalline complexes and densimetry

Dilute aqueous amylose solutions (0.1 wt%) were prepared by dispersing 10 mg amylose in 10 mL ultrapure water. These suspensions were bubbled with nitrogen for 20 min, autoclaved

for 30 min at 160 °C, cooled down to 90 °C and poured into glass tubes containing 0.04 wt% IBU preheated at 90 °C, above its melting temperature (Tm = 76 °C). The tubes were incubated at a fixed temperature (50, 60 or 75 °C). After one week, the tubes were cooled down to room temperature (20 °C). The solid residues were recovered by slow centrifugation (500 rpm, 22 g) and washed three times with water saturated with IBU in order to preserve the crystal stability.

In most cases, recrystallization steps were carried out to obtain larger and more individual single crystals. After the first crystallization, the reaction tubes were incubated at 90 °C for 1 h to partially dissolve the solid fraction, and the recrystallization was carried out for 2-3 days at the temperature used for the first crystallization. This treatment was repeated three times, then the tubes were allowed to cool down. Crystalline complexes were prepared with propan-2-ol for comparison. A DP2500 amylose solution (0.1 wt%) was mixed at 75 °C with 35 vol% preheated propan-2-ol and the solution was kept at 75 °C for one day. The resulting crystals were washed in a 65:35 (vol:vol) water:propan-2-ol mixture [START_REF] Buléon | Single-crystals of V-amylose complexed with isopropanol and acetone[END_REF].

Fragments of a film of VIBU crystals desorbed in a chamber maintaining a 95% relative humidity (r.h.) were floated in 1,2-dichloroethane (d = 1250 g L -1 ) to which 1,1,2-trichloro-1,2,2-trifluoroethane (d = 1560 g L -1 ) was slowly added. When the fragments remained in equilibrium in the mixture, the density was calculated as equal to that of the liquid which was measured using a pycnometer.

X-ray diffraction (XRD)

The crystallized fractions were centrifuged for 10 min at 12000 g. The wet pellets were deposited onto a 7-µm nylon bolting cloth (SaatiTech) and desorbed at 95% r.h. Strips of hydrated mats were introduced into glass capillaries (1.0 or 1.5 mm outer diameter) that were flame-sealed and fixed on collimators (0.25 or 0.50 mm in diameter). The capillaries were X-rayed in a Warhus vacuum camera using a Philips PW3830 generator operating at 30 kV and 20 mA (Ni-filtered CuKα radiation, λ = 0.1542 nm). Two-dimensional diffraction patterns were recorded on Fujifilm imaging plates, read with a Fujifilm BAS 1800-II bioimaging analyzer.

Diffraction profiles were calculated by radially averaging the ring patterns. The sample-to-film distance was calibrated using a calcite standard. The unit cell parameters were calculated with the CelRef module of the LMGP package [START_REF] Laugier | LMGP-Suite: Suite of programs for the interpretation of X-ray experiments[END_REF].

Optical microscopy

Drops of crystal suspensions were deposited on glass slides and covered with glass slips.

The specimens were observed in differential interference contrast mode with a Zeiss Axiophot II microscope equipped with polarized light. Images were recorded with an Olympus SC50 camera.

Transmission electron microscopy (TEM) and electron diffraction (ED)

Drops of dilute crystal suspensions were deposited on glow-discharged carbon-coated grids and allowed to dry. The specimens were observed under low dose illumination with a Thermo Scientific Philips CM200 Cryo microscope operating at 200 kV. For electron diffraction, the TEM grids were mounted on a Gatan 626 holder and fast-frozen into liquid nitrogen just after air drying of the excess liquid. Once in the microscope, the holder was cooled down with liquid nitrogen and the specimen observed at -177 °C. Images and ED patterns were recorded with a TVIPS TemCam F216 camera. ED patterns were recorded from selected areas of about 1 µm 2 and calibrated with a gold-coated carbon film. The intensity of the diffraction spots was measured using the procedure described by Nishiyama et al. [START_REF] Nishiyama | Molecular and crystal structure of 7-fold V-amylose complexed with 2-propanol[END_REF]. In view of the symmetry of the patterns, the intensities of hk0, -hk0, h-k0 and -h-k0 reflections were averaged.

Atomic force microscopy (AFM)

Droplets of dilute crystal suspensions were deposited on pieces of Si wafers and allowed to dry. AFM measurements were carried out with a Molecular Imaging Pico Plus microscope equipped with silicon probes having a nominal spring constant of 37 N m -1 and tip radius of 6 nm (ACT probes from AppNano). The data were recorded in tapping mode in air and analyzed using Gwyddion [START_REF] Nečas | Gwyddion: an open-source software for SPM data analysis[END_REF]. As detailed in Supplementary Data, topographic variations were determined after baseline correction by analyzing the height distribution over different regions of interest (ROIs) in several images (Fig. S2). For a good statistical distribution, a minimum of 20 ROIs were analyzed for each sample. The surface average roughness was calculated over 150 ´ 150 nm 2 regions ( 512´ 512 pixels) of single crystals as:

𝑅 ( = * + , ∑ (𝑍 0 -𝑍 ̅ ) 4 , +
Eq. 1

where 𝑍 ̅ is the average height determined over the ROI.

Attenuated total reflectance Fourier-transform infrared spectroscopy (FTIR)

FTIR spectra of hydrated complexes (equilibrated at 95% r.h.) were recorded with a Perkin-Elmer Spectrum 2 FT-IR spectrophotometer in attenuated total reflection (ATR) mode.

Measurements were done in the 4000-400 cm -1 range with a resolution of 2 cm -1 . Eight repeat scans were averaged for each spectrum.

Solid-state nuclear magnetic resonance (NMR) spectroscopy

Crystal suspensions were centrifuged. The excess liquid was blotted off and the hydrated pellets were packed into a zirconia rotor. Solid-state 13 C NMR analyses were performed with a Bruker Avance III 400 MHz spectrometer ( 13 C frequency of 100.6 MHz) using magic angle spinning (MAS) and cross-polarization (CP). The spinning speed was set at 12 kHz, with a sweep width of 29761 Hz and a recycle delay of 2 s. Each spectrum was averaged over about 6000 scans. The 13 C chemical shifts were calibrated with respect to that of the glycine carboxyl group (176.03 ppm). Single-pulse (SP) solid-state NMR experiments were carried out as well using a 13 C 90° pulse of 4 µs and a delay time of 64 s.

Differential scanning calorimetry (DSC)

Samples (ca. 8 mg) mixed with deionized water (40 µL) were accurately weighed using a Mettler-Toledo XP2U ultra-microbalance and hermetically sealed in large-volume stainless steel DSC pans (Perkin Elmer Pan Kit 0319-1525/1526/1535). Sample pans were scanned against an empty pan as reference in a Discovery differential scanning calorimeter (TA Instruments). The samples were heated from 10 to 150 °C, then cooled down to 5 °C, and immediately reheated to 150 °C, all steps carried out at a rate of 10 °C min -1 .

Molecular modeling and structure refinement

Stereoregular left-handed 7-fold amylose helices with a pitch of 0.80 nm were generated from a α-D-glucosyl residue in the 4 C1 chair shape as described elsewhere [START_REF] Nishiyama | Molecular and crystal structure of 7-fold V-amylose complexed with 2-propanol[END_REF] and energy-minimized using the Universal Force Field [START_REF] Rappé | UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations[END_REF] in the Forcite module of Biovia Materials Studio 5.5 and introduced into the unit cell with parameters and space group determined from the ED and XRD data. IBU and water molecules were distributed in the volumes of the unit cell not occupied by amylose chains and the complete system was energy-optimized again. During this step, the atomic coordinates of amylose were constrained and only IBU and water molecules were allowed to move.

The structure was then refined against experimental ED intensities using the SHELXL program [START_REF] Sheldrick | Crystal structure refinement with SHELXL[END_REF]. The atomic scattering factors for electrons were taken from the International Tables for Crystallography [START_REF] Cowley | Parameterizations of electron atomic scattering factors[END_REF].

Finally, a two-dimensional base-plane ED pattern was calculated from the refined model by convoluting the calculated intensities with a Gaussian function. The similarity between the experimental and calculated patterns was evaluated using the reliability factor R defined as:

𝑅 = 1 -6 𝑰 89: •𝑰 <=>< ‖𝑰 89: ‖ ‖𝑰 <=>< ‖ A 4
Eq. 2

where Iobs and Icalc are the observed and calculated intensities represented as a vector. The components of vector I = (I1, I2, I3, ..., Ii, ..., In) are the intensities of the n reflections considered in the calculation [START_REF] Nishiyama | Molecular and crystal structure of 7-fold V-amylose complexed with 2-propanol[END_REF]. Since the location of the hydrogen atoms could not be determined from the experimental data, hydrogen bonds (H-bonds) were identified from the O...O distance, with a maximum of 0.35 nm. Similarly, the van der Waals interactions were detected from the distance between carbon atoms, with a maximum of 0.40 nm.

Results and discussion

Crystal morphology and influence of amylose DP

In a preliminary series of experiments, the formation of VIBU complexes was investigated for 0.1 wt% amylose solutions containing 0.04 wt% IBU as a function of the 𝐷𝑃 # $$$$$ of amylose and incubation temperature. The DP30 fraction did not form crystalline complexes whatever the temperature. This result agreed with previous observations that a minimum chain length of amylose is required for the formation of stable helical complexes. This value was reported to vary between DP 20 and 40 depending on the complexing agent [START_REF] Gelders | Amylose-lipid complexation: a new fractionation method[END_REF][START_REF] Godet | Crystallization of amylose-fatty acid complexes prepared with different amylose chain lengths[END_REF]. The longer chains (DP200, DP600 and DP2500) formed complexes at all tested temperatures up to 75 °C, while DP60 and DP130

amylose only formed crystals at 50 and 60 °C.

Interference contrast optical micrographs showed that the crystals developed in the form of sea urchin-like spheroidal aggregates of crystals radially growing from a common origin (Figs. 1a andS2). In low-magnification TEM images, slender crystals were seen radiating from thicker electron-opaque regions (Fig. 1b), while flower-like groups of several micrometer-long single lamellae with a more or less well-defined rectangular shape were observed at higher magnification. A higher temperature and shorter amylose chains were preferred to obtain larger and less densely associated single crystals (Fig. S3). However, an additional recrystallization treatment was used to produce more individual single crystals leading to an easier recording of selected area ED patterns. The morphology of the crystals varied with the DP of amylose: the lamellae were longer with the longer amylose chains (Fig. 2a-d) while they were shorter but with a better defined rectangular shape with the shorter chains (Figs. 2e-h and S3). AFM images confirmed the lamellar morphology (Fig. 3a,b) and height distributions were determined over several ROIs (Fig. S4). While the surface roughness was similar for all samples (Rg = 0.7 nm), the lamellar thickness increased with increasing amylose DP (Fig. 3c), reaching a plateau at 7.0 ± 0.4 nm for DP2500, in agreement with the chain folding mechanism described for flexible synthetic polymers [START_REF] Hu | The physics of polymer chain-folding[END_REF] and amylose [START_REF] Manley | Chain folding in amylose crystals[END_REF]Yamashita,1965).

Crystal structure

Base-plane ED patterns were recorded at low temperature from individual frozen-hydrated crystals. The typical diagram shown in Fig. 4 contains about 160 diffraction spots, up to a resolution of 0.34 nm, distributed along a rectangular network defined by two orthogonal reciprocal base vectors a* and b* of very close length. After calibration, the unit cell parameters were found to be a = 2.83 ± 0.01 nm and b = 2.96 ± 0.01 nm. Despite some differences in shape and size, the crystals prepared from all amylose fractions yielded similar ED patterns (Fig. S5a,b). The shape of the crystals, ED diagrams and unit cell parameters are very similar to those of Vpropan-2-ol crystals (Fig. S6a-c) [START_REF] Buléon | Single-crystals of V-amylose complexed with isopropanol and acetone[END_REF][START_REF] Nishiyama | Molecular and crystal structure of 7-fold V-amylose complexed with 2-propanol[END_REF][START_REF] Yamashita | Single crystals of amylose V complexes. II. Crystals with 71 helical configuration[END_REF] and consistent with earlier reports on V-amylose crystallized with different ligands such as thymol, carvacrol [START_REF] Helbert | Données sur la structure du grain d'amidon et des produits de recristallisation de l'amylose[END_REF], geraniol, fenchone, menthone [START_REF] Nuessli | Crystal structure of amylose complexes with small ligands[END_REF], linalool, terpineol [START_REF] Putaux | Single crystals of V-amylose inclusion complexes[END_REF], or pentadecylphenol [START_REF] Kumar | Deciphering structures of inclusion complexes of amylose with natural phenolic amphiphiles[END_REF].

In a typical XRD powder pattern (Fig. 4c) and corresponding profile (Fig. 4d) of hydrated VIBU crystals, the strongest peaks are located at diffraction angles 2θ = 6.9, 8.7, 15.2, 17.7, 18.1, 19.1 and 19.8° (Table S1). The XRD profile is nearly identical to that recorded from the Vpropan-2-ol complex (Fig. S6d) with minor variations in relative peak intensity. After calibration and assuming an orthorhombic structure with P212121 space group, the unit cell parameters were determined as a = 2.824 ± 0.001 nm, b = 2.966 ± 0.001 nm, and c = 0.800 ± 0.001 nm.

No significant difference in terms of crystallinity could be observed between the XRD profiles of complexes prepared with the different amylose fractions (Fig. S5c), in agreement with what was concluded from ED patterns. Moreover, only minor differences in relative intensity of some reflections were observed between the ED patterns (Figs. 4b andS6c) and XRD profiles (Figs. 4d andS6d) recorded from VIBU and Vpropan-2-ol crystals.

FTIR spectroscopy

The FTIR spectrum of VIBU crystals prepared from DP2500 amylose is shown in Fig. 5d and compared to those of pure IBU, pure amylose and a physical mixture of amylose with 10 wt% IBU (Fig. 5a, 5b and 5c, respectively). The carbonyl group (C=O) stretching vibration band of pure IBU observed at 1707 cm -1 is present in the spectra of the physical mixture and VIBU but a new band appeared at 1730 cm -1 in the spectrum of the complex. By analogy with the results from several studies of complexes between IBU and β-cyclodextrins, the band at 1707 cm -1 likely corresponds to a strong interaction between molecules in IBU crystals through H-bonds between carbonyl groups, while the band at 1730 cm -1 would be attributed to less structured H-bonded carbonyls of IBU molecules well-dispersed within the amylose scaffold [START_REF] Braga | Solid state inclusion compound of S-ibuprofen in β-cyclodextrin: structure and characterization[END_REF][START_REF] Bratu | FT-IR spectroscopy on inclusion complexes of β-cyclodextrin with fentbufen and ibuprofen[END_REF][START_REF] Mura | Interactions of ketoprofen and ibuprofen with β-cyclodextrins in solution and in the solid state[END_REF]. The band assigned to aromatic C-C stretching shifted from 1507 cm -1 in pure IBU or the physical mixture to 1512 cm -1 in the complexes. Several bands attributed to IBU that still appeared in the physical mixture, disappeared from the VIBU spectrum. This is the case for absorbances at 669 and 779 cm -1 , attributed to the deformation of the out-of-plane aromatic C-H bonds of the IBU benzene group [START_REF] Ghorab | Elucidation of solution state complexation in wetgranulated oven-dried ibuprofen and β-cyclodextrin: FT-IR and 1 H-NMR studies[END_REF] and bands that correspond to aromatic C-C stretching (1451 cm -1 ), O-H bending (1418 and 934 cm -1 ), C-O stretching (1268-1123 cm -1 ) and in-plane C-H bending (1091-1007 cm -1 ) [START_REF] Coates | Interpretation of infrared spectra, a practical approach[END_REF].

13 C CP/MAS and SP/MAS NMR spectroscopy

The 13 C CP/MAS NMR spectra of hydrated VIBU and Vpropan-2-ol crystals are shown in Figs. 6a andS7, respectively. The peaks in the carbohydrate region (55-100 ppm) can be assigned to the carbons of the glucosyl ring and it was suggested that the position of the carbon C1 resonance was a signature of the number of glucosyl units per turn in the amylose helix (Gidley & Bociek, 1988;[START_REF] Bail | Structural investigation of amylose complexes with small ligands: helical conformation, crystalline structure and thermostability[END_REF]. This resonance located at 103.3 and 103.4 ppm for Vpropan-2-ol and Vmenthone complexes, respectively [START_REF] Bail | Structural investigation of amylose complexes with small ligands: helical conformation, crystalline structure and thermostability[END_REF], was proposed to correspond to a packing of 6-fold helices [START_REF] Biais | Structural and stoichiometric studies of complexes between aroma compounds and amylose. Polymorphic transitions and quantification in amorphous and crystalline areas[END_REF]). However, it was later demonstrated that the Vpropan-2-ol complex actually contained 7-fold helices [START_REF] Nishiyama | Molecular and crystal structure of 7-fold V-amylose complexed with 2-propanol[END_REF]. The C1 resonance of our VIBU and Vpropan-2-ol samples was located at 103. Carignani, Borsacchi, & Geppi, 2011), is absent from the VIBU spectrum, confirming that the excess IBU was properly rinsed out during sample preparation. Instead, the small peak of carbonyl carbons is shifted upfield (173-177 ppm), which is consistent with what was reported for crystalline complexes of IBU and b-CDs [START_REF] Braga | Solid state inclusion compound of S-ibuprofen in β-cyclodextrin: structure and characterization[END_REF]. The region corresponding to aromatic carbons contains a triplet at 141.6, 139.3 and 136.7 ppm, and a doublet at 129.3 and 126.7 ppm. The C6 resonance that characterizes the orientation distribution of the hydroxymethyl group in amylose is maximum at 60.0 ppm for VIBU (Fig. 6) and 60.1 ppm for Vpropan-2-ol (Fig. S7).

Similar values were reported for complexes with propan-2-ol, menthone and linalool that have a similar crystal structure [START_REF] Bail | Structural investigation of amylose complexes with small ligands: helical conformation, crystalline structure and thermostability[END_REF]; Rondeau-Mouro, [START_REF] Rondeau-Mouro | Structural investigation of amylose complexes with small ligands: inter-or intra-helical associations[END_REF].

The SP/MAS NMR spectrum of VIBU was recorded to quantify the guest included in the complex (Fig. 6b). The contribution of IBU is more pronounced than that in the CP/MAS NMR spectrum (Fig. 6a). Such a difference between the CP and SP spectra likely results from a significant mobility of the IBU molecules within the crystal lattice. By deconvoluting the entire C1 resonance between 93 and 110 ppm to separate the contributions from complexed and amorphous amylose [START_REF] Tan | A method for estimating the nature and relative proportions of amorphous, single, and double-helical components in starch granules by C-13 CP/MAS NMR[END_REF], we found that 80 % of the C1 resonance could be attributed to the complex (Fig. S8). Assuming that all IBUs were included in the V-type crystals, the stoichiometry between IBU and glucosyl residues, calculated as the ratio of the integral area of the C1 contribution of amylose (93-110 ppm) to those from aromatic carbons of IBU (123-143 ppm), was calculated to be about 1:8.4. The ratio was 1:7.6

when the region of aliphatic carbons from IBU (16-50 ppm) was considered instead.

Thermal properties

DSC thermograms of VIBU complexes crystallized from a 0.1 wt% DP2500 amylose solution at 75 °C (CIBU = 0.04 wt%) revealed an endotherm with a peak maximum at 130.8 °C with a transition enthalpy of 24.15 J g -1 (Fig. 7). This transition enthalpy is nearly identical to that observed for complexes of a-L-lysophosphatidylcholine with potato amylose leached from common maize starch at saturation [START_REF] Creek | Potential sources of error in the calorimetric evaluation of amylose content of starches[END_REF], suggesting a near maximum complexation with ibuprofen. An exotherm with an onset temperature of 94.2 °C and peak maximum at 83.8 °C was observed on cooling at 10 °C min -1 , with a transition enthalpy of 15.05 J g -1 . Reheating revealed an endotherm at 126.1 °C with a transition enthalpy of 19.4 J g -1 . The rapid crystallization in the calorimeter on cooling resulted in 80% of the complexation yielded by the initial slow crystallization from a dilute solution, and produced a slightly less stable crystalline form with a lower peak endotherm temperature, as could be expected. A shoulder could be seen on both endotherms at about 115 °C, indicating two different thermal events of varying stability. However, the endotherm observed upon reheating the sample after the rapid recrystallization in the DSC was nearly similar to the one recorded during the first heating of the initial complex that was prepared from a dilute solution and with a much slower crystallization. This suggests that the two populations of VIBU complexes do not correspond to different crystal types but rather to IBU molecules with differing interactions with amylose, as one might expect from intra-and interhelical guests. By deconvoluting the endotherm on reheating between 80 and 156.5 °C, it was determined that 87 % of the enthalpy could be attributed to the more stable (i.e. higher temperature) complexed amylose helices, presumably those with intrahelical IBU molecules.

Molecular modeling, structure refinement and interaction analysis

Considering the similarities in morphology, diffraction and spectroscopy features of VIBU and Vpropan-2-ol crystals and since the position of the C1 peak in the 13 C CP/MAS spectrum was consistent with a 7-fold amylose helix, we used the model proposed by Nishiyama et al. for Vpropan-2-ol complexes as a suitable basis to build a plausible model of VIBU. In this orthorhombic unit cell, the helix arrangement can be described by two alternating groups of four left-handed single helices, organized into nearly close-packed and tetragonal motifs, respectively (Fig. S9) [START_REF] Nishiyama | Molecular and crystal structure of 7-fold V-amylose complexed with 2-propanol[END_REF]. This specific organization of 7-fold helices was recently referred to as V7II, as opposed to V7I that would correspond to a pseudo-hexagonal arrangement of closepacked helices [START_REF] Le | Single crystals of V-amylose complexed with bicyclic organic compounds[END_REF]. In Vpropan-2-ol, water and propan-2-ol molecules are located inside the helices and in-between (Fig. S9).

As detailed in the Experimental Section, we have first built a model containing stereoregular left-handed 7-fold amylose helices located at the same positions as those in the Vpropan-2-ol unit cell (Fig. S10a). This amylose scaffold contains two types of voids which dimensions in cross section are suitable to host guest molecules. One corresponds to the columnar cavity inside the helices and the other to the two interhelical channels, both types of cavities being parallel to the c-axis of the unit cell. In this preliminary modeling study, we have only considered the case of (S)-(+)-ibuprofen (S-IBU) as a guest since it is the pharmacologically active form (Scheme S1b). Since one S-IBU molecule is longer than the helical repeat of the amylose helix, a supercell larger than the unit cell along c was built, that contained 3 helical turns, each cavity hosting 2 S-IBU molecules. In each model, consistent with the shape and size of the intra-and interhelical cavities, the guests were broadly aligned with the helical axis of amylose. Inside all helices, the carboxyl group of S-IBU pointed in the same direction as the reducing end of the amylose chain while, in one interhelix cavity, IBUs were all up and in the other, all down.

The relative amounts of S-IBU and water in the supercell were adjusted from the experimental density of the crystals (1.41 g cm -3 measured by densimetry). Since several IBU:water ratios satisfied the density value, we chose the following one: 2 S-IBUs were hosted in each helix while 2 S-IBUs were located in each interhelical cavity. The remaining free interhelical space was filled with water molecules. The supercell thus contained a total of 84 glucosyl residues in the amylose helices, 8 intrahelical and 4 interhelical S-IBUs, and 56 interhelical water molecules, corresponding to a calculated density of 1.41 g cm -3 and an overall stoichiometry of 1 S-IBU for 7 glucosyl units. Although this ratio is a little higher than that determined from the analysis of 13 C SP/MAS NMR data, it corresponds to the maximal loading of the drug in the structure, our goal being to provide a model that can be used later to understand the release profiles of the drug in various media.

The model was first energy-minimized (Fig. S10b) then refined against ED data to a resolution of 0.31 nm, using SHELXL. The final model is shown in Fig. 8 and the corresponding calculated base-plane ED pattern is in very good agreement with the experimental pattern (R = 0.053 calculated with Eq. 2) (Fig. 9). In this structure, the amylose helices have lost their strict stereoregularity as suggested by the dispersion of the torsion angles around the glycosidic linkages and the hydroxymethyl conformation (Table S2), which results

from their local conformational adaptation to the presence of the guest (Table S3). S-IBU molecules deviate from their intrinsically most stable forms established by density functional theory (DFT) modeling and their solution conformation studied by NMR [START_REF] Di Pietro | Assessing the stable conformations of ibuprofen in solution by means of Residual Dipolar Couplings[END_REF]. The conformation of S-IBU can be defined by three torsion angles (Scheme S2). In particular, the T1 and T2 angles describe the orientation of the two substituents with respect to the aromatic ring. In the minima of the internal energy of IBU, the bulky groups carried by the carbons directly linked to the aromatic ring must be distant from the ring, with T1 and T2 close to 120 and 60°, respectively. In our model, they deviate by 30° from their ideal value. Such deviation has also been observed in the crystalline structure of S-IBU [START_REF] Freer | Structure of (S)-(+)ibuprofen[END_REF].

The VS-IBU crystal is strongly stabilized by a large number of intermolecular direct and indirect H-bonds between amylose helices, the latter being mediated by water molecules (Fig. S11a). In addition, the interaction between amylose chains and S-IBU mainly involves dispersion forces, as suggested by the large number of van der Waals contacts (Fig. S11b).

Conclusion

Highly crystalline lamellar complexes were successfully prepared from dilute amylose solutions incubated in the presence of racemic ibuprofen. The single crystal "quality" was optimized by adjusting the amount of added IBU, the incubation temperature, and by performing recrystallization steps. The DP of the amylose fraction (native or biosynthesized in vitro) had an impact on the morphology of the lamellar crystals (size, thickness, aspect ratio) but no significant effect was noted regarding the crystallinity. The base-plane ED and powder XRD data, used to determine the P212121 space group and unit cell parameters of the orthorhombic crystal, along with the 13 C CP/MAS NMR spectrum consistent with the occurrence of a 7-fold amylose helix, suggested that the VIBU and Vpropan-2-ol complexes were isomorphous. Consequently, we proposed a plausible molecular model based on the features of the Vpropan-2-ol complex and that contains four 7-fold antiparallel left-handed single helices per unit cell. IBU molecules are hosted inside the helical cavity but also in the interhelical space along with a number of water molecules. After refinement of the model against ED data, the calculated base-plane ED pattern showed a very good agreement with the experimental pattern, thus validating the 2D projection of the proposed structure, and particularly the distribution of left-handed 7-fold helices within the unit cell. However, the three-dimensional molecular description should still be considered with caution since i) only S-IBU was tested as a guest, ii) no diffraction data along the c-axis was used in the ED-based SHELXL refinement, iii) the stoichiometry of IBU and water in the complex is not known with precision, and iv) for simplicity, the ratio between the number of IBUs and helical repeats was chosen in order to minimize the size of the tested supercell. Considering the assumptions made to build our preliminary model, a more thorough molecular modeling study will be carried out to account for other possibilities regarding the chirality (S or R enantiomer, or both) and the orientation (up or down, both or random) of the guest, as well as the location and stoichiometry of IBU and water molecules.

So far, the structures of the V-type allomorphs that have been solved from diffraction data were always based on left-handed amylose single helices (Brisson, Chanzy, & Winter, 1990[START_REF] Nishiyama | Molecular and crystal structure of 7-fold V-amylose complexed with 2-propanol[END_REF], which raises the question of a possible chiral selectivity of ligands by the amylose scaffold. Indeed, while the chiral selectivity of S-IBU by β-cyclodextrins was demonstrated (Núñez-Agüero, Escobar-Llanos, Díaz, Jaime, & Garduño-Juárez, 2006) and studied for short amylodextrins [START_REF] Hong | Complexation between amylodextrin oligomers and selected pharmaceuticals measured through capillary electrophoresis[END_REF] and maltodextrins [START_REF] Garnero | Ibuprofen-maltodextrin interaction: Study of enantiomeric recognition and complex characterization[END_REF], the question is still open for single-helical amylose in the form of crystalline complexes. Finally, the V7II structural model proposed in this work for VIBU complexes, and in particular, the location of the guest within the amylose scaffold, should serve as a basis to understand the in vitro release properties of ibuprofen in controlled media that mimic the gastro-intestinal environment.

Preparation of amylose by in vitro biosynthesis with amylosucrase -Fractionation and chain length determination:

The reaction was performed at 30 °C in a 50 mM Tris-HCl buffer at pH 8.0 containing 0.5 U mL -1 of purified amylosucrase from Neisseria polysaccharea and 100 mM sucrose Table S3. Conformation of the S-IBU molecules in the asymmetric unit of the VS-IBU complex (two intra-and one interhelical guests). The T1, T2 and T3 torsion angles are defined by the C1-C2-C4-C5, C6-C7-C10-C11, and C7-C10-C11-C12 sequences of atoms, respectively. [START_REF] Freer | Structure of (S)-(+)ibuprofen[END_REF] b conformation corresponding to the energy minimum in the DFT calculation [START_REF] Di Pietro | Assessing the stable conformations of ibuprofen in solution by means of Residual Dipolar Couplings[END_REF]

Figure 1 .

 1 Figure 1. a) Interference contrast optical micrograph and b) low-magnification TEM image of aggregates of VIBU complexes (DP2500 amylose, CIBU = 0.04 wt%, incubation at 75 °C).

Figure 2 .

 2 Figure 2. TEM images of lamellar single crystals of V-amylose complexed with racemic IBU: a,b) DP2500 (native amylose); c,d) DP200; e,f) DP130; g,h) DP60. The incubation temperature was 75 °C for DP2500 and DP200, and 60 °C for the other fractions (CIBU = 0.04 wt%). For all samples, successive recrystallization steps were carried out to obtain more individual crystals.

Figure 3 .

 3 Figure 3. Typical AFM images of lamellar single crystals of DP200 (a) and DP60 (b) amylose complexed with IBU (CIBU = 0.04 wt%, incubation at 75 °C for DP200 and 60°C for DP60); c) average lamellar thickness calculated from height profiles over several regions of interest in the images as a function of amylose 𝐷𝑃 # $$$$$ . The error bars correspond to standard deviations.

Figure 4 .

 4 Figure 4. a) TEM image of a lamellar single crystal of DP60 amylose complexed with IBU; b) corresponding base-plane ED pattern recorded from a frozen-hydrated crystal; c) XRD powder pattern of hydrated VIBU crystalline complexes (DP2500); d) corresponding XRD profile.

Figure 5 .

 5 Figure 5. ATR-FTIR absorbance spectra of racemic IBU (a), native amylose (b), physical mixture of amylose and 10 wt% IBU (c) and VIBU crystalline complexes (d). The series on the right focuses on the 1780-1580 cm -1 region corresponding to the contribution of carbonyl groups.

Figure 6 .

 6 Figure 6. 13 C CP/MAS (a) and SP/MAS (b) solid-state NMR spectra of hydrated VIBU crystalline complexes (DP2500 amylose, CIBU = 0.04 wt%, incubation at 75 °C).

Figure 7 .

 7 Figure 7. DSC thermograms of VIBU complexes crystallized from a 0.1 wt% DP2500 amylose solution at 75 °C (CIBU = 0.04 wt%). The specimen was reheated (heating 2) after the first heating / cooling cycle.

Figure 8 .

 8 Figure 8. a) Axial view of the tentative molecular model of the VS-IBU complex. The position of the 7-fold amylose helices in the unit cell corresponds to that determined for the Vpropan-2-ol complex (Nishyama et al., 2010). S-IBU and water molecules are distributed inside the helices and in-between. The hydrogen atoms have been omitted for clarity. b,c) Longitudinal views of the model highlighting the location and conformation of S-IBU in the intra-and interhelical cavities (b and c, respectively). Amylose helices are drawn as a blue wireframe. The water molecules and hydrogen atoms have been omitted.

Figure 9 .

 9 Figure 9. Comparison of the experimental base-plane ED pattern recorded at low temperature from one VIBU crystal (exp) with that calculated from the VS-IBU model in Fig. 8 (calc).

(

  [START_REF] Potocki-Veronese | Amylose synthesized in vitro by amylosucrase: Morphology, structure, and properties[END_REF]. After 24 h, amylose was precipitated from the reaction medium with 1/4 volume of propan-2-ol at 4 °C for 12 h, and collected by centrifugation (2500 g, 20 min). Part of the total fraction containing chains with a DP ranging from 15 to 220 (𝐷𝑃 # $$$$$ 86)was dissolved in DMSO 10 % (w/v) and was further fractionated by means of preparative sizeexclusion chromatography (SEC) using a 3L column of Biogel P6DG (Biorad). Sugars were eluted with water at 60 °C using a flow rate of 10 mL min -1 . High-performance anion-exchange chromatography with pulsed amperometric detection (HPAEC-PAD) was used to determine the DP distribution of each fraction after dissolution in 1 N NaOH.

Figure S1 .

 S1 Figure S1. HPAEC-PAD chromatograms of the total fraction of amylose (DP86) synthesized in vitro by amylosucrase and the fractions (DP60 and DP130) obtained by fractionation of DP86.

Figure S2 .

 S2 Figure S2. Influence of the DP of amylose on the morphology of the aggregates formed by crystallization with ibuprofen: a) DP2500 (native amylose); b) DP600; c) DP200; d) DP60 (0.1 wt% amylose, CIBU = 0.04 wt%, incubation at 60 °C). Interference contrast optical micrographs of suspensions observed between cover and slip. Scale bar: 20 µm.

Figure S3 .

 S3 Figure S3. Influence of the DP of amylose and incubation temperature on the morphology of the aggregates of VIBU crystals (0.1 wt% amylose solutions, CIBU = 0.04 wt%, first crystallization). TEM images with 1 µm scale bars. "n.c.": no crystallization.

Figure

  Figure S5. a,b) Base-plane ED patterns of VIBU single crystals prepared with DP60 (a) and DP2500 (b) amylose fractions, recorded at low temperature from frozen-hydrated preparations; c) XRD powder profiles of VIBU crystals prepared with various amylose fractions.

Figure S6 .

 S6 Figure S6. a,b) TEM images of lamellar single crystals prepared by crystallizing native amylose in the presence of propan-2-ol; c) base-plane electron diffraction pattern of one crystal recorded at room temperature; d) powder X-ray diffraction profile of a mat of hydrated Vpropan-2-ol crystals.

Figure S7 .

 S7 Figure S7. 13 C CP/MAS NMR spectra of Vpropan-2-ol crystalline complexes (DP2500 amylose).

Figure S8 .

 S8 Figure S8. Deconvolution of carbon C1 peak in the 13 C SP/MAS NMR spectrum of VIBU crystalline complexes (full spectrum in Fig.6). The main peak (blue line) was fitted with a pseudo-Voigt function (purple line) and the amorphous contribution by two Gaussian functions (green and red lines)[START_REF] Tan | A method for estimating the nature and relative proportions of amorphous, single, and double-helical components in starch granules by C-13 CP/MAS NMR[END_REF].

Figure S9 .

 S9 Figure S9. Axial view of the molecular model of the Vpropan-2-ol crystal drawn from the data published by Nishiyama et al. (2010). The hydrogen atoms have been omitted for clarity.

Figure S10 .

 S10 Figure S10. Optimization and refinement steps to produce a tentative molecular model of the VS-IBU complex: a) initial model containing only symmetrical 7-fold amylose helices. The position of the 7-fold amylose helices corresponds to that determined for the Vpropan-2-ol complex (Nishiyama et al., 2010); b) relaxation after addition of S-IBU and water molecules located in the amylose helices and in-between. The conformation of the helices was constrained; c) refinement against ED intensities using SHELXL. All molecules were free to relax. The hydrogen atoms in the displayed models have been omitted for clarity.

Figure S11 .

 S11 Figure S11. Axial views of the VS-IBU molecular model showing the hydrogen bonds (pink dashed lines in a) and van der Waals contacts (blue dashed lines in b). The hydrogen atoms have been omitted for clarity.

Table 1 .

 1 Propan

	-2-ol and racemic ibuprofen

Table 1 .

 1 Molecular characteristics of the amylose fractions used in this study. 𝐷𝑃 # $$$$$ is the weight-average degree of polymerization. 𝑀 # $$$$$ and 𝑀 & $$$$ are the weight-and number-average molecular weights, respectively.

	Fraction	DP30 a	DP60 b	DP130 b DP200 a DP600 a DP2500 c
	𝐷𝑃 # $$$$$$	28	60	130	192	601	2514
	DP min	-	32	32	-	-	-
	DP max	-	123	219	-	-	-
	𝑀 # $$$$$ (g mol -1 )	4500 d	9738	21078	31200 d	97500 d	407268
	𝑀 # $$$$$ 𝑀 & $$$$	1.06 d	1.05	1.07	1.01 d	1.04 d	2.70

).

Table S1 .

 S1 Observed (obs) and calculated (cal) diffraction angles 2θ in the X-ray powder diffraction pattern recorded on hydrated VIBU complexes, and corresponding Miller indices based on an orthorhombic unit cell with a = 2.824 ± 0.001 nm, b = 2.966 ± 0.001 nm, and c = 0.800 ± 0.001 nm.

	h k l	2θ obs (°)	2θ cal (°)
	1 1 0	4.325	4.320
	2 0 0	6.256	6.259
	1 2 0	6.732	6.732
	2 1 0	6.934	6.933
	2 2 0	8.656	8.646
	1 3 0	9.467	9.477
	3 1 0	9.856	9.857
	2 3 0	10.93	10.923
	3 2 0	11.110	11.132
	1 4 0	12.336	12.341
	2 0 1	12.724	12.724
	4 1 0	12.894	12.888
	3 4 0	15.208	15.210
	5 2 0	16.831	16.795
	3 5 0	17.685	17.671
	5 3 0	18.080	18.087
	6 1 0	19.056	19.090
	6 2 0	19.785	19.788
	5 5 0	21.731	21.725
	3 7 0	23.022	23.017
	4 7 0	24.485	24.500
	7 4 0	25.131	25.127

Acknowledgements

We thank LabEx ARCANE and CBH-EUR-GS (grant agreement #ANR-17-EURE-0003) for funding the PhD thesis of C.A.K.L, the Glyco@Alps program (grant agreement #ANR-15-IDEX-02), the NanoBio-ICMG Platform (FR 2607, Grenoble) for granting access to the AFM, NMR and TEM facilities, and Henri Chanzy (CERMAV) for stimulating discussions and critical reading of the manuscript.

Yang, L., Zhang, B,. Yi, J., Liang, J., Liu, Y., & Zhang, L.-M. (2013). Preparation, characterization, and properties of amylose-ibuprofen inclusion complexes.

Starch/ Stärke, 65, 593-602. Zhang, L., Cheng, H., Zheng, C., Dong, F., Man, S., Dai, Y., & Yu, P. (2016). Structural and release properties of amylose inclusion complexes with ibuprofen. Journal of Drug Delivery Science and Technology, 31,[101][102][103][104][105][106][107]

Supplementary Data

Scheme S1. Chemical structure of amylose (a), (S)-(+)-ibuprofen (b) and (R)-(-)-ibuprofen (c), drawn with MarvinSketch (https://chemaxon.com/products/marvin).

Determination of the average lamellar thickness and rugosity from AFM images:

First, large ( 35´ 35 µm 2 ) and high-resolution ( 4096´ 4096 pixels) images of the sample were recorded (Fig. S2a). Smaller regions of interest ( 4´ 4 µm 2 , 470 ´ 470 pixels, without large aggregates) were selected from these images for height distribution analysis (Fig. S2b).

Special care was taken for background plane fitting of each region. The peaks in the height distributions were fitted with Gaussian functions, allowing to determine the height of the Si support, the individual lamellae, as well as stacks of two or more crystals (Fig. S2c). The height of one lamella was calculated by averaging the height of individual and stacked lamellae after subtraction of the height of the support.


Sabine Montaut
email: smontaut@laurentian.ca

Benjamin S Guido

Claude Grison

Patrick Rollin

Mey

Identification of Glucosinolates in Seeds of Three Brassicaceae Species Known to Hyperaccumulate Heavy Metals

Keywords: glucosinolate, Alyssum fallacinum, Noccaea caerulescens, Iberis intermedia, LC-MS

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

Introduction. -The long history of mining operations has led to the accumulation of trace elements (TE) in the environment. TE are persistent in ecosystems and living organisms.

As they are not biodegradable, they tend to concentrate easily in living organisms along food chains in the magnification process. The high toxicity of TE in soil, water resources, and crops affects public health. In spite of their toxicity, heavy metals can also exert a selective pressure on living organisms and thus drive evolution. Metal-tolerant plant species are able to grow on metal-contaminated soils while metal-hyperaccumulating plant species can extract, transport, and concentrate metals from soils into their above-ground parts.

About 450 metalhyperaccumulators have been discovered throughout the world. Among them, metalhyperaccumulating plants belonging to the Brassicaceae family have been the most used to develop wide programs for the phytoremediation of contaminated sites. Recently, they have gained considerable interest because of their potential recovery in green chemistry. Metal-rich biomass allows the production of new catalysts, referred to as ecocatalysts. Ecocatalysts provide increased yields in chemical production and increased regio-and chemoselectivity, which results in high added value [1] [2]. This new approach to using metal-rich biomass, such as Alyssum spp. and Thlaspi (Noccaea) spp. could spur the development of phytoextraction, a technique considered promising for long, yet without viable economic outlets [3].

Production of a group of plant secondary metabolites called glucosinolates (GLs) is a common feature to the Brassicaceae family [START_REF] Blažević | Glucosinolates[END_REF]. Discussions on the role of GLs in defence and tolerance mechanisms in metal hyperaccumulators are often based on total GL content.

However, the relationship between GLs and metal accumulation, and these relationships regarding defence against herbivores in hyperaccumulator species, are not clear. Toward a better understanding of these relationships, we chose to investigate the GL profile of seeds of three Brassicaceae species known to hyperaccumulate heavy metals, using LC-MS. We investigated A. fallacinum is a known hyperaccumulator of Ni [START_REF] Cecchi | [END_REF][6][7]. I. intermedia, commonly known as a variety of candytuft, is a known hyperaccumulator of Tl [8][START_REF] Al-Najar | and sustainability of agro-ecosystems through basic and applies research[END_REF][START_REF] Scheckel | [END_REF]. The GLs of I. intermedia and A. fallacinum have never been investigated. N. caerulescens is a known hyperaccumulator of Zn, Cd, and Ni [11][12][13][14][15][16][17][18]. This plant is also Pb-hypertolerant but not necessarily Pbhyperaccumulating [19]. Glucoiberin (1), glucosinalbin (6), gluconapin, 4hydroxyglucobrassicin (7), and 4-methoxyglucobrassicin were identified in N. caerulescens shoots (5 weeks old) by liquid chromatography-atmospheric pressure chemical ionization mass spectrometry of the desulfo-GLs [13] [20]. Gluconasturtiin can also be present in shoots [13].

Furthermore, gluconapin, 6 (major), and sometimes neoglucobrassicin were found in the roots.

In this study, glucoputranjivin, glucomalcolmiin, 1, and glucocapparin were also detected in leaf tissue [11]. In addition, sinigrin (3), 6, and 4--rhamnosyloxybenzyl GL (8), also known as glucomoringin, were extracted from seeds and leaves [21].

The selection of the hyperaccumulator plants mentioned above was undertaken in the context of our ongoing phytoremediation programs in Southern Europe [3].

<insert the chemical formula about here> fallacinum seed extract displayed two major peaks at t R 6.4 and 20.9 min (Fig. 1a). The compound at t R 6.4 min had identical t R , mass (422 a.m.u), and UV spectra as a commercial standard of 1. The compound at t R 20.9 min, had a mass of 406 a.m.u and a UV spectrum similar to that of a commercial standard of glucoerucin. Therefore, 2 was identified as glucoibervirin.

The LC-MS chromatogram of A. fallacinum from our study contained the same major GLs 1 and 2 as those found in the seeds of Alyssum peltarioides Boiss. and Alyssum sibiricum Willd. from Turkey [22].

Glucosinolate Composition of Iberis intermedia seed. Four GL peaks were observed in the chromatogram of the I. intermedia seed extract (Fig. 1b). 1 and 2 were also identified in this extract. The peak at t R 8.7 min with a mass of 358 a.m.u. was determined to be 3 by comparison to a commercial standard. The peak at t R 23.4 min with a mass of 408 a.m.u was found to be glucotropaeolin (4). This identification was confirmed by comparison of the UV and mass spectra and t R of a commercial standard of 4. Iberis umbellata L. seeds, purchased commercially, were previously found by GC-MS analysis of GL hydrolysis products to contain 1-3 [22]. The major GL was 1 (8.4 mmol/100 g sample) followed by 3 (3.9 mmol/100 g sample), and 2 (0.14 mmol/100 g sample) [22]. In the leaves, 3 and 4 (83 mmol/10 g dried leaves, 1 mmol/10 g dried leaves, respectively) were found and quantified [23]. Other Iberis spp. are also known to contain 1 and 2 [22,[24][25][START_REF] Kjaer | Fortschritte der Chemie Organischen Naturschtoffe[END_REF][START_REF] Jaki | [END_REF].

Glucosinolate Composition of Noccaea caerulescens seed. The chromatographic profiles of N. caerulescens were quite similar (Fig. 2a and2b). Glucoconringiin was determined to be compound 5 at t R 8.5 min with a mass of 390 a.m.u by comparing UV and mass spectra and t R with an authenticated sample previously isolated in our group from Bretschneidera sinensis Hemsl. seeds [28]. Another minor compound at t R 17.9 min, with a mass of 424 a.m.u was identified as 6 by comparing its UV, mass spectra, and t R with those of a commercial standard.

Compound 7 was eluted at t R 20.8 min and had a mass of 463 a.m.u. The spectroscopic data were similar to those of 7 previously found in our group in Brassica elongata Ehrh. seeds [29].

In addition, the major peak at t R 21.3 min with a mass of 570 a.m.u was identified as 8 by comparison of its UV, mass spectra and t R with an authenticated sample [30]. All other unidentified peaks in all extracts possessed UV spectral characteristics of flavonoids [START_REF] Markham | The flavonoids[END_REF].

We report for the first time the presence of 5 in N. caerulescens. However, 5 was previously deduced from the detection of 5,5-dimethyloxazolidine-2-thione, the myrosinase hydrolysis product of 5 in Thlaspi kovatsii Heuff. (synonym Thlaspi avalanum Pančić) and Thlaspi alpestre Jacq. (synonym Noccaea alpestris (Jacq.) Kerguélen) seeds [22]. Contrary to previous studies in seeds [21], we did not detect 3 which was claimed the major GL in some accessions. In our case, the major GL was 8, which was also the case in other accessions from Spain (Valle de Varrados and Pontaut), Luxembourg (Lellingen), and France (Navacelles) [21].

Our results confirm intraspecific variation in the GL profile in N. caerulescens seeds. In addition, 7, which had only been detected in the shoots [20,21], was detected in the seeds for the first time. We did not detect 1, 4-methylsulfanylbutyl GL, glucoputranjivin, glucomalcomiin, nor glucocapparin, which were previously reported in leaf tissue [11].

A previous investigation has shown that when N. caerulescens was exposed to increased Zn concentrations, the concentration of 6 in shoots diminished, whereas in the roots, the concentration of 6 increased with Zn concentration [13]. Another study confirmed that GL concentration in shoots decreased when the concentration of Zn in leaves increased [16].

However, a separate research group demonstrated that GL production increased in shoot tissue with Zn treatments, especially 4-methylsulfanylbutyl GL, also known as glucoerucin [11]. Other researchers indicated that the total GL content in Noccaea leaves increased with the concentration of Ni or Cd in the soil [12]. However, the total GL amount in damaged leaves was higher than in undamaged leaves in the presence of high concentrations of Ni [12]. Conversely, it was shown that the GL content was higher in undamaged leaves in the presence of a high concentration of Cd [12]. Another investigation reported a higher GL content in undamaged than in damaged leaves of N. caerulescens grazed by thrips (Frankliniella occidentalis) and grown in Zn-contaminated soil while the anthocyanin concentration was high in damaged leaves [START_REF] Asad | [END_REF]. Furthermore, GLs were shown to deter gastropods from eating N. caerulescens [15] [18].

Moreover, the damage caused by the cabbage whitefly (Aleyrodes proletella) did not affect the production of GLs in N. caerulescens grown in Zn-contaminated soil [16]. In addition, N.

caerulescens was shown to contain organic acids, amino acids, metallothioneins, and phytochelatins which could be responsible for its heavy metal-hyperaccumulating capacity [11] [14]. Finally, galactolipids, anthocyanins, nicotianamine, and oxylipins were also found in the plant [11] [32].

Conclusions. -Our investigation of the three plants known to hyperaccumulate heavy metals A. fallacinum, I. intermedia, and N. caerulescens by LC-MS has shown the efficacy of the method in separating their GLs. We report for the first time the presence of 1-2 in A. fallacinum and 1-4 in I. intermedia. Finally, the GL profile including 5-8 in N. caerulescens from two sites of collection showed no differences but some discrepancies with literature data were pointed out. H O to 100% MeOH over 60 min [37]. After each run, the initial mobile phase conditions were and the system was allowed to equilibrate. The flow rate was kept constant at 1 mL min -1 .

The column temperature was held at room temperature. The HPLC was interfaced to an Agilent

  Alyssum fallacinum auct. non Hausskn. (synonym Alyssum baldaccii Vierh. ex Nyár.), Iberis intermedia Guers., and Noccaea caerulescens (J. Presl & C. Presl) F.K. Mey. (synonym Thlaspi caerulescens J. Presl & C. Presl).

  Results andDiscussion. -The seeds of three plants known to hyperaccumulate heavy metals (A. fallacinum, I. intermedia, and N. caerulescens) of the Brassicaceae family were extracted and analysed by LC-MS for intact GLs. We collected samples from one site each for the first two plant species from sites known to be contaminated with heavy metals (see Experimental Part). N. caerulescens seeds were collected from two different contaminated sites. Glucosinolate Composition of Alyssum fallacinum seed. The chromatogram of the A.

Fig. 1 .Fig. 1 .

 11 Fig. 1.

Fig. 2 .

 2 Fig. 2.

Fig. 2 .

 2 Fig. 2. HPLC chromatograms of the methanolic extract of Noccaea caerulescens seeds from a)

  Glucoerucin and glucoiberin were purchased from Cfm Oskar Tropitzsch (Marktredwitz, Germany). Glucotropaeolin was purchased from Chromadex (Irvine, CA, USA). Sinigrin was purchased from Sigma Aldrich (Oakville, ON, Canada). HPLC-grade H 2 O was generated in the laboratory through a Nanopure Diamond Ultrapure water system by Barnstead (Dubuque, IA, USA).HPLC-ESI-MS Analysis.The extracts were dissolved in MeOH/H 2 O 7/3 (v/v) (2.5 mL for A. fallacinum, 6 mL for I. intermedia, and 5 mL for N. caerulescens from Bergenbach and Les Avinières) and were filtered through a plug of cotton prior to analysis by a high-performance liquid chromatograph (HPLC). The analyses were performed by injecting 10 µL of extract into an Agilent Technologies HP 1100 (New Castle, DE) HPLC equipped with a quaternary automatic injector, diode-array detector (wavelength range 190-600 nm), degasser, and a Hypersil ODS column (5 m, 4.6  200 mm). The two mobile phase solvents, MeOH and H 2 O, were prepared with 0.15% triethylamine and 0.18% formic acid added as ion-pairing reagents.Both solutions were filtered using 0.45 m nylon membranes. The initial mobile phase was 100% HPLC-grade H 2 O. At 10 min, the mobile phase was switched to a linear gradient of 100%

	Plant Material. A. fallacinum seeds were harvested in 2014 in Anogia-Gonies road, 1 km
	from Sisorha towards Gonies (estimated 35°17.9' N, 24°55.7' E (Greece)), on a serpentine soil
	containing high concentrations of Ni (1,350 ppm), and identified by Prof. Alan J.M. Baker
	(School of Botany, University of Melbourne, Australia). I. intermedia seeds were harvested in
	2010 in Les Avinières, Saint-Laurent-Le-Minier, Gard (03°66'50''E, 43°93'13'' N, France)
	which is a mining site in which the soil contains Zn (up to 156,000 ppm), Pb (36,354 ppm), Cd
	(700 ppm), and Tl (115.1 ppm) [33] [34]. N. caerulescens seeds were harvested in 2012 in Les
	Avinières (same mining site where I. intermedia seeds where collected) and 2010 in Bergenbach,
	sges (France) on a serpentine soil containing Ni (116 mg kg -1 ammonium acetate-EDTA
	ractable element), Zn (25 mg kg -1 ammonium acetate-EDTA extractable element) and Cd
	0 mg kg -1 ammonium acetate-EDTA extractable element) [35] [36]. I. intermedia and N.
	caerulescens were identified by Prof. Claude Grison (University of Montpellier, France).
	Extract Preparation. Seeds (361 mg of A. fallacinum, 512 mg of I. intermedia, 575 mg
	of N. caerulescens from Bergenbach, and 517 mg of N. caerulescens from Les Avinières) were
	Experimental Part frozen in liquid N 2 and ground to powder with a mortar and pestle. The powder was extracted
	General. All solvents were ACS grade and used as such. Formic acid was purchased with boiling EtOH/H 2 O (7/3 v/v) (2  5 mL) for 5 min. The solutions were concentrated to
	from BDH (Toronto, ON, Canada). HPLC-grade MeOH, absolute EtOH, and triethylamine dryness (38 mg of A. fallacinum, 89 mg of I. intermedia, 62 mg of N. caerulescens from

The financial support from the National Sciences and Engineering Research Council of Canada (NSERC Research Tool and Instruments), Canadian Foundation for Innovation (Leaders Opportunity Fund)-Ontario Research Fund, Laurentian University, Employment and Social Development Canada (Canada Summer Jobs) as well as from the Centre National de la Recherche Scientifique (CNRS, France) is gratefully acknowledged. The authors also thank Prof. A.J.M. Baker for discussions and identification of A. fallacinum. (reagent grade) were purchased from Fisher Scientific (Whitby, ON, Canada). Glucosinalbin was purchased from Apin Chemicals Ltd. (Abingdon, UK). Bergenbach, and 70 mg of N. caerulescens from Les Avinières).


Benjamin Deneu

Alexis Joly

Pierre Bonnet

Maximilien Servajean

François Munoz

How do Deep Convolutional SDM trained on satellite images unravel vegetation ecology?

Keywords: Species Distribution Model, Convolutional Neural Network, Ecological Interpretation, Plant functional traits, Trait-environment relationships, t-SNE

Species distribution models (SDM) assess and predict how species spatial distributions depend on the environment, due to species ecological preferences. These models are used in many different scenarios such as conservation plans or monitoring of invasive species. The choice of a model and of environmental data have strong impact on the model's ability to capture important ecological information. Specifically, state-of-the-art models generally rely on local, punctual environmental information, and do not take into account environmental variation in surrounding landscape. Here we use a convolutional neural network model to analyze and predict species distributions depending on high resolution data including remote sensing images, land cover and altitude. We show that the model unravel the functional response of vegetation to both local and large-scale environmental variation. To demonstrate the ecological significance of the results, we propose an original statistical analysis of t-SNE nonlinear dimension reduction. We illustrate and test the traitsspecies-environment relationships learned by the model and expressed in t-SNE dimensions.

Introduction

Understanding and predicting the geographic distribution of species is a key objective in ecology and conservation. Species Distribution Models (SDM) aim at characterizing the relationship between the environment and species occurrences, depending on species ecological niches [START_REF] Elith | Species Distribution Models: Ecological Explanation and Prediction Across Space and Time[END_REF]. The ecological niche is multidimensional, and involves factors playing in a complex fashion (i.e. non linear system) and at multiple spatial scales. Most SDMs are correlative methods relating known species occurrence data to potential environmental predictors [START_REF] Guisan | Predictive habitat distribution models in ecology[END_REF][START_REF] Antoine | Predicting species distribution: offering more than simple habitat models[END_REF][START_REF] Peterson | Ecological niches and geographic distributions[END_REF].Popular examples of such methods include MAXENT [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF][START_REF] Phillips | Modeling of species distributions with maxent: new extensions and a comprehensive evaluation[END_REF], random forest [START_REF] Cutler | Random forests for classification in ecology[END_REF] and boosted regression trees [10,[START_REF] Elith | A working guide to boosted regression trees[END_REF]. Earlier works devised SDMs based on single-layer neural networks [START_REF] Lek | Application of neural networks to modelling nonlinear relationships in ecology[END_REF][START_REF] Baran | Stochastic models that predict trout population density or biomass on a mesohabitat scale[END_REF], and recently deep neural networks have proved relevant to better address the complexity of ecological niches, allowing to recognized a larger complexity in the way environment shapes ecological niches [START_REF] Chen | Deep multi-species embedding[END_REF][START_REF] Benkendorf | Effects of sample size and network depth on a deep learning approach to species distribution modeling[END_REF]. Key advantages of deep learning are that (i) it allows characterizing complex structuring of ecological niche depending on multiple environmental factors, (ii) it can learn ecological features common to a large number of species, and thus grasp the signatures of common ecological processes and improve SDM predictions across species [START_REF] Botella | A deep learning approach to species distribution modelling[END_REF][START_REF] Christin | Applications for deep learning in ecology[END_REF].

A specific class of neural networks initially proposed in [START_REF] Lecun | Backpropagation applied to handwritten zip code recognition[END_REF], named Convolutional Neural Networks (CNN), has very recently been proposed for SDM [START_REF] Botella | A deep learning approach to species distribution modelling[END_REF][START_REF] Deneu | Location-based species recommendation using co-occurrences and environment-geolifeclef 2018 challenge[END_REF]. A property of CNN is that they rely on spatial environmental tensors rather than on local values of environmental factors. These tensors represent the spatial realization of environmental factors around each point. Unlike other SDM approaches, CNNbased SDMs (CNN-SDMs) can use this very large input data and therefore potentially capture richer information than in punctual vectors. CNNs were originally designed for image classification [START_REF] Lecun | Backpropagation applied to handwritten zip code recognition[END_REF] and proved to outperform any other statistical or machine learning methods in the task of learning complex visual patterns. This largely explains why CNN-based architectures are the most popular deep learning implementations in ecological studies since few years [START_REF] Christin | Applications for deep learning in ecology[END_REF].

In this paper, we use a deep species distribution model based on a convolutional neural network applied to high resolution remote sensing images as one of the input covariables. This model has already been evaluated and compared to other more classical methods revealing superior performance to state of the art models. In this paper we are interested in the interpretation of these performances. The major contribution of the paper is to provide an ecological interpretation of the model. Therefore, we propose to rely on t-SNE [START_REF] Maaten | Visualizing data using t-sne[END_REF], a nonlinear dimension reduction technique widely used in data science to visualize high-dimensional feature spaces. More precisely, we use it to build a low-dimensional embedding of the 2028-dimensional representation space learned by the convolutional neural network (i.e. the 2028 neurons used as input of the species prediction layer of the model). The resulting lowdimensional embedding is then used in two ways: (i) to visualize the relationships between the learned features, the species traits and the environment, and (ii) to quantify these relationships using statistical tests. The main outcome of our study is that the statistical tests clearly demonstrate that the model was able to capture meaningful relationships between the species traits and the environment. This is particularly remarkable in the sense that none of these information were used as input covariables during the training of the model. The model was actually capable of inferring them directly from the high-resolution spatialized images used as input and the species occurrences.

Material and Methods

CNN-SDM Model Training and Validation

Training Dataset A detailed description of the dataset (named GeoLifeCLEF 2020) is provided in [START_REF] Cole | The geolifeclef 2020 dataset[END_REF]. It consists of 1,921,123 observations from the US (1,097,640) and France (823,483) covering 31, 435 plant and animal species (Table 1). Each species observation is paired with high-resolution covariates (RGB-IR imagery, land cover and altitude), see Fig. 1. We used RGB-IR imagery from the 2009-2011 cycle of the National Agriculture Imagery Program (NAIP) in the U.S. 6 , and from the BD-ORTHO® 2.0 and ORTHO-HR® 1.0 databases from IGN in France 7 . Land cover data originates from the National Land Cover Database (NLCD) [START_REF] Homer | Completion of the 2011 national land cover database for the conterminous united states-representing a decade of land cover change information[END_REF] for the U.S. and from CESBIO 8for France. All elevation data comes from the NASA Shuttle Radar Topography Mission (SRTM) 9 . All of these high-resolution covariates were homogenized at a spatial resolution of 1 meter per pixel and provided as 256 × 256 images covering a 256m × 256m square centered on each observation (some were oversampled and some where downsampled). Deep Convolutional SDM Architecture Our deep convolutional neural network is the composition of non linear transformations of the input space z = φ(x) with a linear classifier ψ(z). The vector z is called the representation vector, or feature vector of the input tensor x. The architecture used is based on the Inception V3 [START_REF] Szegedy | Rethinking the inception architecture for computer vision[END_REF] model adapted in the same way introduced in [START_REF] Deneu | Participation of lirmm/inria to the geo-lifeclef 2020 challenge[END_REF]. The representation layer has size 2048 and the loss is the cross-entropy. The output of the model can be interpreted as the probability for each species.

Predictive Performance We evaluated the model and compared it to other more classical approaches, such as a random forest learned on environmental data (RF), using a spatial block holdout procedure 10 . On a top-30 score, the CNN obtained 23.5% success compared to the 20.4% success of the random forest environmental model. This evaluation allowed us to validate the performance of the model (CNN-SDM) that we will study later in this report.

Ecological Interpretation of the Learned Features

We considered environmental and species trait data not used during model training. We assess the extent to which the model could capture ecological information such as the functional response of plant species to environmental constraints (specifically, climatic and pedological factors). In the following subsections, we first describe the environmental and species traits data used (2.3). In subsection 2.4, we then describe the non linear dimension reduction technique that we used to embed the 2048-dimensional feature vectors z into a low-dimensional space of 2 or 3 dimensions aimed at conserving only the most structuring information. In subsection 2.5, we describe the statistical tests that were performed on top of the resulting lowdimensional embedding for the ecological interpretation. Finally, in subsection 2.6, we describe some additional methodological details that were used for visualization purposes.

Environmental and Trait Data

We used 19 bio-climatic rasters (30arcsec 2 /pixel (above 1km 2 /pixel) from WorldClim [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF]) and 8 pedologic rasters (250m 2 /pixel, from SoilGrids [START_REF] Hengl | Soilgrids250m: Global gridded soil information based on machine learning[END_REF]), the detailed list and resolutions are presented in Table 2.

We used ecological information available for more than 1,400 plant species, in terms of Ellenberg's indicator values (EIVs) [START_REF] Ellenberg | Vegetation ecology of central Europe[END_REF] (see Table 3). These variables represent an ordinal classification of ecological strategies for major environmental constraints and the use of essential resources [START_REF] Bartelheimer | Functional characterizations of ellenberg indicator valuesa review on ecophysiological determinants[END_REF]. 10 Test occurrences are contained in 5×5km quadrats with no train occurrences and represent 2.5% of the overall set.

Dimension Reduction

CNNs train a representation space in which species occurrences projections z = φ(x) tend to be linearly separable (thanks to the multinomial logistic regression ψ trained on top of z). We analyzed how the structure of the learned space could grasp major ecological and environmental information unused during learning. We identified major dimensions of the learned space by projecting the feature vectors z into a very low-dimensional space of 2 or 3 dimensions. For this purpose, we used the t-SNE dimension reduction method [START_REF] Maaten | Visualizing data using t-sne[END_REF]. First, t-SNE constructs a probability distribution over pairs of high-dimensional objects such as similar objects are assigned a higher probability while dissimilar points are assigned a lower probability. Second, t-SNE defines a similar distribution over the points in the low-dimensional map, and it minimizes the Kullback-Leibler divergence (KL divergence) between the two distributions with respect to the locations of the points in the map. The main advantage of t-SNE is that it is able to preserve the main regularities of the original space even if some of them are not linear. We also computed PCA as a preliminary dimension reduction step.

In more details, to process the dimension reduction we randomly selected 32, 000 training occurrences x j and computed their representations z j = φ(x j ). Then, we first reduced the dimension from 2048 to 50 by PCA using the PCA method of the scikitlearn package. The resulting 50-dimensional feature vectors were then further reduced by t-SNE using the t-SNE method of the scikitlean package. For most experiments, the used dimension for t-SNE was set to 2 (apart from the geographical map of Figure Figure 3b where it was set to 3 and base on other occurrences, see section subsection 2.6). In the following, we denote as z = g (z) the resulting 2dimensional feature vectors where the function g denotes the complete dimension reduction function (PCA+t-SNE).

Ecological Interpretation of t-SNE Dimensions

Relationship of t-SNE Dimensions to Ecological Traits and Environmental Factors

For each of the two variables z1 and z2 , corresponding in the two axis of the t-SNE, we fitted a linear model using either the environmental variables or the species trait variables as input covariates (lm function in R environment). The link between the representation space learned by the model (condensed in the t-SNE space) and these data (not used during the learning process) offers a way to interpret the ecological significance of patterns learned by the CNN-SDM models.

Traits-Environment Relationships

Trait-environment relationship represent essential functional responses of plants to changing environmental conditions. Characterizing these relationships is an essential goal in functional biogeography [START_REF] Violle | The emergence and promise of functional biogeography[END_REF]. We tested whether t-SNE axes unravel such relationships by mean of weighted correlation analysis. We calculated the correlation between species trait and environmental value for each training occurrence, weighted by the score of the occurrence along each t-SNE axis. Under the null hypothesis, the t-SNE axis does not reflect any linkage between trait and environment, under the alternative hypothesis, increasing score is associated to joint variation of trait and environment. We measured the weighted correlation [START_REF] Miller | Functional traits and community composition: a comparison among community-weighted means, weighted correlations, and multilevel models[END_REF] for each t-SNE variable zk and for each pair of ecological trait and environmental variable. We used the function wtd.cor package weights in R. Furthermore, because we computed many tests for all pairs of trait and environmental variables, we corrected the p-value of the tests by method false discovery rate (p.adjust function in R).

Visualization

We used different types of visualisations to illustrate the main sources of ecological variation expressed in the two-dimensional t-SNE space. We discretized the space of the t-SNE by applying an n ×n-sized grid mesh. For each cell of this grid we retrieved the m-coordinate vector of the center of the cell and we search for the nearest neighbor among the set of vectors z using scikitlearn's nearest neigbhor function. Once the closest point to the center of each cell is retrieved, only those that are actually inside the respective cell are kept by a filter on these coordinates. This results in associating for each cell the occurrence closest to the center if there is one. The different figures of this visualization consist in the display, for each cell, of a data relative to the associated occurrence. For the first fig. 3a, we displayed the RGB patch of the remote sensing imagery corresponding to this occurrence. For the other figures (2a and 2b) we displayed the value of an ecological trait of the corresponding species or the realization of one of the environmental variables at the point of occurrence.

The second representation is based on the 3-dimensional t-SNE space and, conversely, represents the realization of this space on the geographical space. This 3dimensional t-SNE space was obtained by the same process as the 2-dimensional space described in detail in section 2.4 but with one difference, the dimension reduction was not performed with the same set of occurrences. As this representation is geographical, it requires a better coverage of the territory than the simple random selection of a certain number of occurrences. We selected occurrences with a method of nearest neighbor. For a grid of points with a resolution of 1km over the whole territory, we associated each point to the nearest occurrence. The selected occurrences are those used for dimension reduction. Then to plot the 3-dimensional t-SNE on the map, we re-scale each axis such as the values are included in 0-255. Each point can then be associated to a RGB color by its coordinates. The resulting map is the plot of each point as a color pixel at its geographical position.

Results and Discussion

The statistical analysis of the learned t-SNE space reveals that the CNN-based SDM could grasp important ecological processes shaping the large scale distributions of plants. Indeed, almost all linear models linking scores along t-SNE to species traits provided highly significant relationships (Table 4). We can in particular note a strong coefficient on EIV T (species preferences in Temperature) on both axes. The plant Ellenberg indicator values could alone explain a large part of the variation of species scores along t-SNE axes (Adjusted R 2 of 0.111 and 0.231). These results show that the model is able to capture information related to the ecology of the species in the patches. Fig. 2a represent the species preferences in temperature and Fig. 2b represent Annual mean temperature of occurrences data on the t-SNE. We found a clear gradient of mean annual temperature on the t-SNE. The fact that this gradient is particularly sharp on t-SNE, which is the result of a dimension reduction, indicates that information strongly correlated with annual temperature is captured and important in the model. It is quite logical that the figure on the associated ecological trait (EIV T) is similar confirming that the plants with the coolest temperature preferences are found in the coldest areas and vice versa. These results highlight that the information captured by the model is strongly related to the environment and ecology of the species, even though this data was not used directly. Fig. 3a represent RGB patches from the model training data on the t-SNE space and Fig. 3b represents the geographical projection of the 3-dimensional t-SNE space on the territory with the colors defined by the position in the t-SNE space.

These two figures highlight two complementary pieces of information. On the one hand, the display of the patches on the t-SNE highlights the landscape factors that the CNN is able to identify and differentiate in its representation space. Indeed the Fig. 3a highlights different areas of the t-SNE corresponding to several major types of distinct landscapes. We can note for example the mountain areas on the left next to which we find the forests a little lower down. The more agricultural areas also stand out while the lower and right part of the t-SNE is dominated by urban type landscapes. On the other hand, Fig. 3b shows learning on a larger scale. This map shows how the CNN has learned to distinguish large eco-geographic regions. Of particular note is the Mediterranean region (on southern part of the map), and the mountainous areas (i.e. the Pyrenees on the south-western part, and the Alps on the south-eastern part), which seem to be well identified. More generally the map allows to show coherence on a large scale, between typical vegetation zones well identified by ecologists, and bio-geographic zones well characterized on the field in previous studies. The confrontation of these two results shows that the model is capable of discerning both large regions and different habitats within the same large geograph-ical area (such as urban-rural). The ability to identify factors at different scales based on fine scale imaging training data is an important result that highlights one of the advantage of the CNN models. Concerning the species trait-environment relationships, the detailed results matrix of each pair for each axis is given in supplementary materials Tables C1 and C2 [START_REF] Deneu | Supplementary materials: How do deep convolutional sdm trained on satellite images unravel vegetation ecology?[END_REF]. A general finding is that most of the species traitenvironment correlations weighted by the t-SNE axes are highly significant. This shows that the CNN is capable of learning a relationship between environment and species ecology through the use of high-resolution covariate patches. All these results highlight a strong point of the CNN model which is capable of capturing information related to the environment and the ecology of species through the use of data accessible at large scale and fine resolution, but also to highlight links between the environment and the ecology of species. This model could therefore be useful in many ecological scenarios where the objective would be to study this link. All these points are also insights for interpreting the performance of the CNN models. By learning a representation space common to all the studied species, the model is able to extract common information between species that are directly correlated to environmental factors or related to the ecology of the species.

Conclusion

Our study shows the interest of using convolutional neural networks for species distribution modelling (CNN-SDM), at high spatial resolution over large geographical areas. The use of the t-SNE technique made it possible for the first time to visualise and test the learning capacities of this type of model to capture relationships between ecological strategies of plants and environmental conditions at their occurrence locations. The methodological framework established here offers a new way of statistically assessing the extent to which the relationships between the plant traits and their environments are correlated, inspired by the conceptual and methodological framework of functional biogeography [START_REF] Violle | The emergence and promise of functional biogeography[END_REF]. It also allows visualizing through t-SNE the clusters of visual information (from satellite images) deemed most relevant during the learning phase of the CNN-SDMs. This work shows that the approach captures richer information on the landscape context than the specific ponctual data related to the environment for predicting the presence of species, and will undoubtedly open up new perspectives in the analysis of plant-environment trait relationships.

  Example of high-resolution covariates patch (respectively: RGB imagery, IR imagery, altitude, land cover).

Fig. 1 :

 1 Fig. 1: Occurrences distribution (training data in blue, test data in red) in France (a), and an example high-resolution tensor (b). RGB patch is displayed with native colors, IR imagery, altitude and land cover are in artificial colors from purple (lowest values) to yellow (highest values).

Fig. 2 :Fig. 3 :

 23 Fig. 2: Visualization of the occurrences' species preferences in temperature (EIV T)(a) and the annual mean temperature (bio_1) at occurrences location (b) on the t-SNE space. Artificial colors from purple (lowest values) to yellow (highest values).

Table 1 :

 1 Number of occurrences, of species and of occurrences per species grouped for plants and animals, and per region (for the training set).

	Nb

. of occurrences Nb. of species Nb. of occ. per sp.

	Kingdom	US	France	US	France	US	France
	Plants	524,280	741,010	11,369	3,114	46.1	238.0
	Animals	551,563	61,865	13,882	4,899	39.7	12.6
	All	1,622,120	802,875	25,251	8,013	42.6	100.2

Table 2 :

 2 Summary of environmental variable rasters provided.

	Name Description

Table 3 :

 3 Summary of the used plant species traits.

	Name	Description	Ranges of values (Nb species)
	EIV L	Light availability	2 -9 (1,423)
	EIV T	Temperature	1 -9 (1,413)
	EIV K	Climatic continentality	1 -8 (1,411))
	EIV AirH Air humidity	1 -9 (1,405)
	EIV F	Soil humidity	1 -12 (1,405)
	EIV R	Reaction (soil acidity / pH) 1 -9 (1,410)
	EIV TroL Trophic level	1 -9 (1,412)
	EIV S	Salt (soil salinity)	0 -9 (1,416)
	EIV SoiT Soil texture	1 -9 (1,416)
	EIV N	Nitrogen (soil fertility)	1 -9 (1,422)

Table 4 :

 4 Ellenberg linear models

		tsne_1	tsne_2
		Estimate (Std. Error)	Estimate (Std. Error)
	EIV L	-1.562 * * * (0.236)	2.630 * * * (0.214)
	EIV T	-5.035 * * * (0.260)	9.478 * * * (0.235)
	EIV K	2.523 * * * (0.389)	0.150 (0.352)
	EIV AirH	0.732 * (0.378)	-0.869 * * (0.341)
	EIV F	1.461 * * * (0.409)	2.328 * * * (0.370)
	EIV R	0.461 * * (0.213)	-0.686 * * * (0.192)
	TroL	2.902 * * * (0.181)	3.851 * * * (0.163)
	EIV S	-2.714 * * * (0.366)	-0.981 * * * (0.331)
	EIV SoiT	1.649 * * * (0.252)	2.848 * * * (0.228)
	EIV N	0.211 (0.212)	-0.598 * * * (0.192)
	Constant	-5.951 (3.947)	-113.831 * * * (3.569)
	R 2	0.111	0.232
	Adjusted R 2	0.111	0.231
	Note:

* p<0.1; * * p<0.05; * * * p<0.01

National Agriculture Image Program, https://www.fsa.usda.gov

https://geoservices.ign.fr

http://osr-cesbio.ups-tlse.fr/~oso/posts/2017-03-30-carte-s2-2016/

https://lpdaac.usgs.gov/products/srtmgl1v003/

Acknowledgement

This project has received funding from the French National Research Agency under the Investments for the Future Program, referred as ANR-16-CONV-0004 and from the European Union's Horizon 2020 research and innovation program under grant agreement No 863463 (Cos4Cloud project).