Dr Carole Mathe
email: carole.mathe@univ-avignon.fr

Elodie Mezzatesta

Nathalie Dupuy
email: nathalie.dupuy@univ-amu.fr

Francesca Modugno
email: francesca.modugno@unipi.it

Ludovic Bellot-Gurlet
email: ludovic.bellot-gurlet@sorbonne-universite.fr

Salvado Nativitat
email: nativitat.salvado@upc.edu

Keywords: SIMPLISMA, MIR, PCA, Ancient-Egypt, Mummification, Balms SIMPLISMA, MIR, PCA, Ancient-Egypt, Mummification, Balms. 26 27

This study presents the efficiency of an original method to characterize mummification balms. Currently, long protocols are used to extract each substance from mummification balms by using techniques such as gas chromatography and mass spectrometry. The SIMPLISMA approach using a simple algorithm combined with chemometric treatment, principal component analysis, was used to interpret the 47 Fourier Transform Infrared spectra. This treatment gave 6 pure variables and their corresponding extracted concentration profiles. The first pure spectrum corresponded to residual species, the second allowed the identification of polycyclic aromatic hydrocarbons, molecules which have a pyrolytic origin. Proteins were characterized, and they could be considered as markers of a good conservation of individuals. Finally, two pure spectra showed the same natural substance, fat, but on the one hand unaltered and on the other hand degraded fatty matter. This fast, simple and noninvasive approach allowed the composition of each balms studied as well as their state of conservation/degradation to be characterized.

Introduction 28

Balms in the mummification process are used for thousands of years in the history of ancient 29

Egypt. Given the importance of this funerary practice, the study and chemical characterization 30 of balms take on their full meaning, but also allow them to be placed in an Egyptian context. [START_REF] Kaiser | An index of factorial simplicity[END_REF] is therefore interesting to study each substance used in the formulation of the balm. A wide 32 variety of natural substances is used in the preparation of balms, such as plant oils, animal 33 fats, resins, beeswax and even other additives such as henna or spices [START_REF] Goyon | Rituels funéraires de l'ancienne Egypte. Le rituel de l'Embaumement. Le 314 rituel de l'Ouverture de la Bouche. Les livres des Respirations[END_REF][START_REF] Nicholson | Ancient Egyptian Materials and Technology[END_REF][START_REF] Buckley | Complex organic chemical balms of Pharaonic 319 animal mummies[END_REF][START_REF] Abdel-Maksoud | A review on the materials used during the 321 mummification processes in ancient Egypt[END_REF][START_REF] Marshall | Les momies egyptiennes. La quête millénaire d'une 324 technique[END_REF]. [START_REF] Magi | Distribution of polycyclic aromatic 407 hydrocarbons in the sediments of the Adriatic Sea[END_REF] substances have certain properties such as being hydrophobic, antibacterial or odoriferous . 35 Each material used in the formulation of a balm is subject to a succession of chemical reactions 36 linked to its evolution over time, from the moment it is obtained to the ageing of the balm, but 37 it can also undergo reactions linked to the preparation, such as heating. Compounds of 38 alteration can have an important part in the identification of plant species. The molecular 39 composition of samples can therefore provide information on the nature of the constituents, 40 the state of conservation of the sample and the treatment undergone by the material during 41 the preparation and/or the preparation of balms [START_REF] Peters | [END_REF][START_REF] Evershed | Organic residue analysis in archaeology: the archaeological biomarker 329 revolution[END_REF][START_REF] Ménager | Analyse de baumes issus de tête momifiée 331 (Thèbes)[END_REF]. 42

Different analytical techniques are used for the chemical characterization of balms depending 43 on the nature of the constituents. Inorganic materials are analyzed by using spectroscopic 44 tools such as Raman spectroscopy, energy dispersive X-ray spectroscopy and X-ray 45 diffraction. Gas chromatography coupled to mass spectrometry is mainly used for the analysis 46 of the organic part. The use of this separative technique is often preceded by a pre-treatment 47 of the samples, such as saponification reactions, extractions with hexane, dichloromethane or 48 diethyl ether and derivatization reactions (e.g. trimethylsilylation) [START_REF] Colombini | Characterization of the Balm of an 333 Egyptian Mummy from the Seventh Century B.C[END_REF][START_REF] Łucejko | An 336 analytical approach based on X-ray diffraction, Fourier transform infrared spectroscopy 337 and gas chromatography/mass spectrometry to characterize Egyptian embalming 338 materials[END_REF][START_REF] Ménager | Study of egyptian mummification balms by 340 FT-IR spectroscopy[END_REF][START_REF] Jacqueline | Les produits d'embaumement égyptiens: nouvelles données 343 pharmacologiques[END_REF][START_REF] Goyon | A propos de la momie infantile du musée 345 Calvet d'Avignon[END_REF][START_REF] Brettell | Organic residue 347 analysis of Egyptian votive mummies and their research potential[END_REF][START_REF] Łucejko | Chemical analyses of 350 Egyptian mummification balms and organic residues from storage jars dated from the 351 Old Kingdom to the Copto-Byzantine period[END_REF][START_REF] Jones | A prehistoric Egyptian mummy: Evidence for an 'embalming 355 recipe' and the evolution of early formative funerary treatments[END_REF]. Chromatographic 49 analyzes often provide a fine chemical characterization but it is also a laborious and 50 discriminating technique. Fourier transform infrared spectroscopy (FT-IR) is also used in this 51 archaeological context. It is a fast and relevant vibrational technique for the analysis of both 52 organic and inorganic matters, i.e. it allows the mineral part to be characterized and different 53 classes of organic materials to be highlighted within the same sample [START_REF] Łucejko | An 336 analytical approach based on X-ray diffraction, Fourier transform infrared spectroscopy 337 and gas chromatography/mass spectrometry to characterize Egyptian embalming 338 materials[END_REF][START_REF] Daher | Advanced discriminating 358 criteria for natural organic substances of Cultural Heritage interest: Spectral 359 decomposition and multivariate analyses of FT-Raman and FT-IR signatures[END_REF]. FT-IR requires 54 moreover a small amount of sample to obtain an imprint of the overall chemical composition 55 of each individual balm and FT-IR can be non-destructive and non-invasive [START_REF] Daher | Towards a non-invasive quantitative analysis of 362 the organic components in museum objects varnishes by vibrational spectroscopies: 363 Methodological approach[END_REF]. This 56 spectroscopic tool is increasingly used for the study of cultural heritage artefacts because of 57 all of these advantages and thus it is a technique of choice for analyzing mummy balms. 58

This global analytical technique is nevertheless restrictive due to the limited identification of 59 the main components. The obtained spectra are the result of the additivity of all the spectral 60 bands of all the molecules present. Owing to the complex chemical composition of 61 mummification balms, in this study a curve resolution method for self-modeling analysis named 62 SIMPLISMA is used. This mathematical method allows pure spectra and associated 63 "contributions" of different substances present in balms without prior information of their 64 existence to be obtained. 65

This study deals with the ATR-FTIR spectra obtained of 47 balms samples from 35 human 66 Egyptian mummies belonging to the collection of the Confluence Museum (Lyon, France). The 67 aim of this study was to determine a fine chemical composition of the balms in terms of fatty 68 matter, beeswax, bitumen, polysaccharides, proteins, resinous material and their degradation 69 Confluences Museum in Lyon (France) were considered (Table 1). The mummies studied were 78 female and male dating from Predynastic period to the Roman Period and coming from 79 different necropolises such as Thebes, Esna, Kom-Ombo and Deir-el-Medineh. The weight of 80 the samples varied from 20 to 200 mg and they were collected at different location points on 81 the mummies: on the mummified skin, on the textile covered by balms, or directly in the balms 82 mass. 83 84 2.2. Infrared spectroscopy 85 FT-IR analyzes were performed on a Nicolet™ iS10 with a MCT detector cooled with liquid 86 nitrogen. Spectra were acquired in Attenuated Total Reflectance mode with a diamond crystal 87 and using Omnic Nicolet software. 100 scans of symmetrical interferograms were averaged 88 and the spectrum was calculated from 1800 to 700 cm -1 at 4 cm -1 spectral resolution. Each 89 analysis was repeated for take into account the heterogeneity of the samples. 90 91 2.3. Chemometrics 92

SIMPLISMA 93

The method used for self-modeling mixture analysis the SIMPLISMA method described in the 94 literature [START_REF] Gemperline | Mixture analysis using factor analysis I: Calibration and quantitation[END_REF][START_REF] Hamilton | Mixture analysis using factor analysis. II: Self-modeling 367 curve resolution[END_REF][START_REF] Windig | Interactive self-modeling mixture analysis[END_REF][START_REF] Windig | Self-modeling mixture analysis of second-derivative near-372 infrared spectral data using the SIMPLISMA approach[END_REF][START_REF] Windig | The use of second-derivative spectra for pure-variable based self-modeling 375 mixture analysis techniques[END_REF][START_REF] Windig | Spectral data files for self-modeling curve resolution with examples using 378 the Simplisma approach[END_REF]. This interactive method is used for self-modeling mixture analysis by 95 resolving mixture data into pure component spectra and concentration profiles without the help 96 of prior information about the mixture. When overlapping spectral features are present in 97 spectroscopic data, this tool is unable to resolve broad spectral components and separate 98 spectral absorption bands characterizing one component. Its concept is based on the 99 determination of pure variables (e.g., a wavelength or a wavenumber in spectroscopic terms) 100 that have received the contributions from only one component and an optimization using least-101 square method [START_REF] Windig | Simple-to-use interactive self-modeling mixture analysis of FTIR 381 microscopy data[END_REF]. In order to obtain the proper resolution of the mixture data, user interaction 102 is required to deal properly with features such as noise, peak shift, and instrument drift [START_REF] Currie | Metrological measurement accuracy: discussion of "measurement error 384 models[END_REF]. 103 Some applications of SIMPLISMA were already reported in the literature, e.g., on Fourier 104 transform Raman spectra for the time resolved activation of hydrogen peroxide by a nitrile [START_REF] Vacque | Self-Modeling Mixture 386 Analysis Applied to FT-Raman Spectral Data of Hydrogen Peroxide Activation by 387 Nitriles[END_REF], 105 gamma irradiated polymers [START_REF] Gaston | Evaluation of multilayer film stability by 390 Raman spectroscopy after gamma-irradiation sterilization process[END_REF], FT-IR microscopy spectra of a polymer laminate, and 106 pyrolysis mass spectra of biomaterials (feedstock) [START_REF] Gemperline | Mixture analysis using factor analysis I: Calibration and quantitation[END_REF][START_REF] Hamilton | Mixture analysis using factor analysis. II: Self-modeling 367 curve resolution[END_REF][START_REF] Windig | Self-modeling mixture analysis of second-derivative near-372 infrared spectral data using the SIMPLISMA approach[END_REF][START_REF] Windig | Simple-to-use interactive self-modeling mixture analysis of FTIR 381 microscopy data[END_REF][START_REF] Windig | The use of second-derivative spectra for pure-variable based self-modeling 375 mixture analysis techniques[END_REF][START_REF] Windig | Spectral data files for self-modeling curve resolution with examples using 378 the Simplisma approach[END_REF]. The fitting of the 107 SIMPLISMA results is calculated using the Relative Residual Sum of Squares (RRSSQ), 108 square sum between the calculated and the original spectra. The SIMPLISMA data treatments 109 were performed with Matlab 7.8 software. The up-ward shift of the baseline due to the 110 fluorescence is corrected with Savitzky Golay derivative pre-processing. 111

PCA analysis 112

PCA is a descriptive statistical method whose purpose is to present, in graphical form, the 113 maximum amount of information contained in a data table with the values of p quantitative 114 variables for n individuals. In this study, the algorithm used to calculate eigenvalues is NIPALS 115 (Non-linear Iterative Partial Least Squares) developed by Herman Wold in 1966 [START_REF] Wold | Path Models with Latent Variables: The NIPALS Approach[END_REF]. The 116 infrared data were reprocessed by the Unscrambler software by applying certain corrections 117 to eliminate instrumental variations. The first one is a correction of the baseline because the 118 variations of the baseline participate in a variance and therefore, they can affect the different 119 spectra. The second correction is normalization (Standard Normal Variate SNV) which allows 120 the intensity of each spectra is at the maximum to be checked. 121 Kaiser-Meyer-Olkin test is used to validate sampling adequacy [START_REF] Dziuban | When is a correlation matrix appropriate for factor 396 analysis? Some decision rules[END_REF]. Kaiser, H.F, mentioned 122 that "the test measures sampling adequacy for each variable in the model and for the complete 123 model. The statistic is a measure of the proportion of variance among variables that might be 124 common variance" [START_REF] Kaiser | An index of factorial simplicity[END_REF]. The lower the proportion, the more suited the data is to Factor 125 Analysis. KMO returns values between 0 and 1. A rule of thumb for interpreting the statistic 126 KMO values less than 0.6 indicate the sampling is not adequate. In our PCA, the KMO equal 127 0.773 [START_REF] Cerny | A Study Of A Measure Of Sampling Adequacy For Factor-401 Analytic Correlation Matrices[END_REF]. 128 for ten minutes and then centrifuged at 6000 rpm for ten minutes. The supernatant was 139 removed, and this step was repeated twice. The three fractions were combined and evaporated 140 to dryness under a stream of nitrogen. Trimethylsilylation was realized with N,O-141 Bis(trimethylsilyl)trifluoroacetamide 99%/Trimethylchlorosilane 1% (BSTFA/TMCS) at 70°C for 142 thirty minutes. The mixture was evaporated to dryness under a stream of nitrogen, it was then 143 dissolved in 0.2 to 1 mL of hexane/dichloromethane (1/1, v/v) and finally it was filtered on a 144 0.45 µm PTFE cartridge filter (VWR). 1 μL of the solution was injected into GC-MS. 145 146

GC-MS conditions 147

A Thermo Scientific Focus gas chromatographic system composed by a Thermo Scientific Al 148 3000 auto-sampler coupled with an ITQ 700 ion trap mass spectrometer was used (Thermo 149 Fisher Scientific). The GC column was a fused silica capillary column Thermo trace GOLD TG-150 5MS (5% diphenyl / 95% dimethylpolysiloxane, 30 m length × 0.25 mm i.d. × 0.25 μm film 151 thickness). Helium was the carrier gas with a constant flow of 1 mL.min -1 . 1 μL of each sample 152 was injected in splitless mode. The injector temperature was set at 250 °C. Mass spectra were 153 recorded in electron impact mode with an electron ionization voltage of 70 eV, an ionization 154 time of 25,000 μs and a mass range of 40-650 m/z. The transfer line, ion trap and manifold 155 temperatures were respectively set at 300 °C, 200 °C and 50 °C. Detected compounds were 156 identified from interpretation of their mass spectra with reference to NIST database, retention 157 time and published literature. The oven temperature started at 50 °C with a 2 min hold, increase 158 of 10°C.min -1 up to 200 °C, increase of 2.5 °C.min -1 up to 310 °C, increase of 8°C.min -1 up to 159 330 °C and held isothermally for 3 min. The total analysis time is 67 minutes. 160 161

Results and discussion 162

The first part of this work was to determine the characteristic bands of the main chemical 163 functions of the various natural substances which could be present in mummification balms 164 and that could be characterized in the mid-infrared. These substances were identified using 165 literature and fresh reference materials. As a reminder, the analyses were carried out between 166 1800 and 700 cm -1 , this zone was chosen because from it is the fingerprint region. The IR 167 assignments of the bands are provided in supplementary material according to the literature 168 (Socrates, 2001). Despite their different locations on the human mummies and the number of 169 deceased individuals studied, balms samples have similar MIR spectra. A broad band (3700-170 3000 cm -1 ) is due to OH stretching vibrations (νOH) (oleuropein, cellulose, organic acids, etc.) 171 and the two bands (3000-2800 cm -1 ) from symmetric and asymmetric CH stretching vibrations 172 (νCH, alkyl). The 1800-1500 cm -1 region corresponds to C O and C C stretching 173 vibrations (esters, acid, carboxylate, proteins, aromatic ring). The 1500-1200 cm -1 range is 174 very complex with especially CH (δCH, alkyl), and OH deformation vibrations as CO stretching 175 vibrations (phenols). The intense bands between 1150 and 950 cm -1 corresponds mainly to 176 the endocyclic and exocyclic C-O stretching vibrations of carbohydrates. Absorption bands in 177 the range 900-746 cm -1 (C-H) are characteristic of the number of adjacent hydrogen atoms 178 on the aromatic ring. 179 180 Each FT-IR experiment was realized on 47 balms samples in duplicate so the SIMPLISMA 181 treatment was applied to the 94 spectra. This treatment gave 6 pure variables and their 182 corresponding extracted concentration profile. A pure variable is defined as a variable whose 183 intensity mainly results from one of the components in the mixture under consideration [START_REF] Windig | Interactive self-modeling mixture analysis[END_REF]. 184

Extracted contribution profiles do not represent a real concentration but the contribution of 185 each pure component in each spectrum. Concentration profiles at zero do not indicate an 186 absence of this substance but they rather traduce that it is not the primary information within 187 the spectrum of this sample, meaning that this pure spectrum is surely characterized by 188 another substance. This phenomenon can be explained by the fact that the infrared is an 189 additive analysis of the data. 190

The first pure spectrum corresponded to instrumental contribution (baseline, water 191 contribution) as well as residual species such as the mineral, so it is not more detailed in this 192

study. 193

The second pure spectrum (Fig 1 .A) was formed by taking the band at 736 cm -1 as the main 194 one. Other bands of lower intensity could also be observed at respectively 1480, 1420, 931, 195 867 and 809 cm -1 . All of these bands correspond to the vibrations of the C-H aromatic bonds 196 [START_REF] Ménager | Analyse de baumes issus de tête momifiée 331 (Thèbes)[END_REF]. These bonds were characterized by numerous bands depending on the nature of the C-H. 197

In view of the bands present which are characteristic of aromatic C-H bonds and the total 198 absence of bands characteristic of the presence of oxygen atoms, this pure spectrum could be 199 attributed to pure Polycyclic Aromatic Hydrocarbons (PAHs). These compounds can appear 200 after a reaction of carbonization and their presence have been already described in balms [START_REF] Ménager | Analyse de baumes issus de tête momifiée 331 (Thèbes)[END_REF]. 201

The sample concentration profile corresponding to this pure PAHs' spectrum was in large 202 proportion in the samples referred 44, 45, 46, 52 and 63 (Fig 1 .B). These results were 203 confirmed by chromatographic analyzes. PAHs molecules were detected in these samples by 204 gas chromatography coupled to mass spectrometry. In order of elution, PAHs compounds 205 were identified to naphthalene, acenaphthene, phenanthrene, anthracene, fluoranthene, 206 pyrene, triphenylene, benzanthracene, perylene, benzopyrene, benzoperylene and 207 dibenzopyrene, with a relative proportion of 95% in all of the corresponding chromatograms. 208

This relative percentage was calculated from the area of all the peaks detected. The origin of 209 these PAHs can be determined by considering the peak area ratios of phenanthrene 210 (Phe)/anthracene (Ant) and fluoranthene (Fluo)/pyrene (Pyr). These four molecules are stable 211 over time and undergo, in an archaeological context, few variations in their concentrations. 212 Thus, the ratios Phe/Anth > 30 and Fluo/Pyr<1 indicate a petrogenic origin of PAHs and the 213 ratios Phe/Anth <10 and Fluo/Pyr>1 indicate a pyrolytic origin of PAHs. These molecules are 214 notable markers of weathering processes [START_REF] Perra | Polycyclic aromatic hydrocarbons 404 pollution in sediments: distribution and sources in a lagoon system (Orbetello, Central 405 Italy)[END_REF], for example, under high temperature 215 phenanthrene will form anthracene and thus this phenomenon conduces to decrease the ratio 216 Phe/Anth [START_REF] Magi | Distribution of polycyclic aromatic 407 hydrocarbons in the sediments of the Adriatic Sea[END_REF]. Regarding all of the ratios calculated, PAHs come from a combustion of organic 217 matter, therefore a pyrolytic weathering process of anthropogenic origin. However, no visible 218 traces of combustion can be observed on the deceased individuals. Two hypothesis can be 219 formulated: either the deceased individuals were burned before wrapping, or they were already 220 wrapped before burning and wrapped again. 221

The third pure spectrum extracted (Fig 2 .A) showed several intense bands around 1650 and 222 1550 cm -1 corresponding to the vibrations of the amide functions which may come from 223 proteins [START_REF] Barth | Infrared spectroscopy of proteins[END_REF]. The concentration profile associated with this spectrum (Fig 2 .B) showed that 12 224 samples [START_REF] Peters | [END_REF][START_REF] Ménager | Analyse de baumes issus de tête momifiée 331 (Thèbes)[END_REF][START_REF] Ménager | Study of egyptian mummification balms by 340 FT-IR spectroscopy[END_REF][START_REF] Goyon | A propos de la momie infantile du musée 345 Calvet d'Avignon[END_REF][START_REF] Daher | Advanced discriminating 358 criteria for natural organic substances of Cultural Heritage interest: Spectral 359 decomposition and multivariate analyses of FT-Raman and FT-IR signatures[END_REF][START_REF] Hamilton | Mixture analysis using factor analysis. II: Self-modeling 367 curve resolution[END_REF][START_REF] Windig | Interactive self-modeling mixture analysis[END_REF][START_REF] Magi | Distribution of polycyclic aromatic 407 hydrocarbons in the sediments of the Adriatic Sea[END_REF]50,61, 64 and 68) contained a significate proportion of 225 proteins. The presence of proteinous material was confirmed by specific Amido Black 10B tests 226 [START_REF] Mathe | Compréhension des mécanismes de coloration des liants 412 protéiques picturaux à l'aide du Noir Amide 10B[END_REF]. There are different possible sources of protein such as the human body of the deceased, 227 an additional ingredient added in the balm or an external contamination. However, this 228 substance was found in high proportion in the samples that had the least amount of alteration. 229 Proteins could therefore be markers for the well preservation of the mummy. 230

The fourth pure extracted spectrum (Fig 3 .A) showed an extremely intense band at 1050 cm -1 231 which was characteristic of the presence of polysaccharides. Figure 3.B showed the 232 concentration profile of each sample for the pure polysaccharide spectrum. Polysaccharides 233 could be obtained from a variety of sources. They could come from the textile of the strips, the 234 gommo fraction of a resinous material or also from a gum like Arabic gum. The presence of 235 polysaccharides in all samples was therefore perfectly explainable. 236

The two last pure extracted spectra and their two respective concentration profiles were 237 presented in Figure 4. Both of these pure extracted spectra (Fig. 4A and4C) characterized the 238 same natural substance (i.e. fatty matter or wax, such as beeswax) but in different state of 239 conservation. Figure 4A presented a high intensity band at 1704 cm -1 (Fig 4 .A) which was often 240 observed in mummy balms corresponding to the C=O function of the carboxylic acids. This 241 chemical group was generally found in the free fatty acids coming from altered and/or degraded 242 fatty matter. In the figure 4C, the first observable band was at 1740 cm -1 , characteristic of the 243 vibrations of the C=O bond of the esters that could come from the triglycerides of plant oils, 244 animal fats or waxes as well as beeswax [START_REF] Svečnjak | An 414 Approach for Routine Analytical Detection of Beeswax Adulteration Using FTIR-ATR 415 Spectroscopy[END_REF]. A band was present at 1168 cm -1 , interpreted to 245 the C-O bond of the ester and carboxylic acid functions [START_REF] Pelletier | Modification et photopolymérisation d'huiles végétales en vue de leur 417 application dans les encres et vernis UV[END_REF]. Moreover, other bands were 246 observed at respectively 1463, 1378 and 722 cm -1 corresponding to C-H bonds of aliphatic 247 compounds. These bonds come from the long carbon chains of beeswax and plants oil or 248 animal fats in the form of triglycerides [START_REF] Cuní | Characterization of the binding medium 420 used in Roman encaustic paintings on wall and wood[END_REF]. Both concentration profiles presented 249 heterogeneous positive and negative profiles with similar intensities (Fig. 4B and4D). All 250 samples contained fatty matter or beeswax well conserved or not, exception to those 251 containing PAHs or another predominant substance. These results were illustrated by Principal 252

Component Analysis (PCA) and these associated loadings in order to facilitate their readability. 253 PCA indicated that 2 principal components (PCs) explained 79% of the variability with 55% for 254 PC1 and 24% for PC2. The obtained results were illustrated in the figure 5 and figure 6. Three 255 distinct groups have been formed representing samples with a significate proportion of 256 altered/degraded (group No. 1) or well conserved (group No. 2) fatty matter or beeswax and 257 the group no 3 contained samples representative by proteins, polysaccharides… However, 258 exceptions were observed for the samples referred 6, 39, and 50. 259 GC-MS analyses allowed these interpretations to be corroborated. The group No. 1 included 260 samples with a high proportion of short carbon chain monocarboxylic free fatty acids and many 261 dicarboxylic fatty acids (Table 2). The formation of short-chain fatty acids and dicarboxylic 262 acids corresponds to a dehydration reaction of fatty matter. This reaction is natural and 263 corresponds to rancidity [START_REF] Fritsch | A study of volatile compounds produced by the 422 autoxidation of methyl oleate, oleic acid and cis-9-octadecene[END_REF]. It is an oxygenation of unsaturated fatty acids by oxygen from 264 the air. Several radical chemical reactions form a hydroperoxide intermediate to ultimately 265 generate short monocarboxylic acids (C) and dicarboxylic acids (DC). The dicarboxylic acids 266 formed generally have between 4 and 13 carbons; azelaic acid DC 9 is the acid mostly formed. 267

The diacids formed depend on the position of the unsaturation of the fatty acids initially present 268 [START_REF] Passi | Saturated 425 dicarboxylic acids as products of unsaturated fatty acid oxydation[END_REF]. The group No. 2 included samples with few degradation marker, i.e. dicarboxylic fatty 269 acids (Table 3). These groups were therefore a good match for the respective concentration 270 profiles of the degraded (group No. 1) and the undegraded fatty matters (group No. 2). 271 272

Conclusion 273

Owing to the various natural substances present in mummification balms associated to several 274 phenomena of alteration and/or degradation, the interpretation of FT-IR results needs great 275 care. In addition, due to the chemical complexity of these archaeological samples, it can be 276 difficult to characterize all their materials. In this study, the SIMPLISMA algorithm was applied 277 with success to extract and to characterize several ingredients used in balms formulation. It 278 was possible to identify the presence of PAHs, proteins, polysaccharides and fatty matter or 279 beeswax well conserved or degraded/altered. Moreover, it was also possible to consider and 280 to compare a significate number of samples with a few experiments and to compare them. GC-281 MS was performed to corroborate and to validate the previous interpretations. The obtained 282 extracted pure spectra associated to their respective concentration profiles provided precious 283 and fine results. Infrared analyses coupled with SIMPLISMA and PCA analysis also allowed 284 the presence a set of chemical families present in complex matrices such as balms to be 285 determined. This work allowed various substances present in mummies' balms to be 286 highlighted. PAHs, the only compounds that do not correspond to natural substance but they 287 have an anthropic origin. Mummy balms containing PAHs were apart and they were pointed 288 out by SIMPLISMA. However, some substances were not detected at all: mineral matter, 289 bitumen, and resins. Middle FT-IR is not the method of choice for mineral analysis. Bitumen is 290 a non-major constituent in these balms and is already difficult to detect without an extraction 291 protocol specific to this chemical family. Plant resins are however in great quantity in certain 292 balms but this product has no characteristic bands. Resinous materials have common bands 293 with the free fatty acids. For all of these reasons, it was not possible to discriminate them and 294 mathematically the algorithm did not manage to extract the resinous substances. 295

This study presented the efficiency of the method and it proposed an original approach to 296 characterize mummification balms. It could be interesting to use a portable infrared device, 297 which is non-invasive and requires no sampling, as a first approach with the parallel use of 298 SIMPLISMA. This approach would give a clear idea of the composition of the balms of very 299 This research did not receive any specific grant from funding agencies in the public, 310 commercial, or not-for-profit sectors. 311

;0 9:0, 10:0, 14:0, 16:0, 17:0, 18:1 (tr), 18:0, 20:0, 22:0, 24:0 5, 6, 7, 8, 9, 10 26:0, 28:0, 30:0 9 8;0 9:0, 10:0, 12:0, 14:0, 15:0, 16:0, 17:0, 18:0 7, 8, 9 24:0, 26:0, 28:0 10 8:0, 9:0, 10:0, 14:0, 15:0, 16:0, 17:0, 18:1, 18:0 5, 6, 7, 8, 9, 10 -14 8:0, 9:0, 10:0, 14:0, 15;0, 16:0, 18:1, 18:0, ricinoléique 5, 6, 7, 8, 9, 10 -15 8:0, 9:0, 10:0, 12:0, 14:0, 15:0, 16:0, 18:0 7, 8, 9, 10 -22 8;0 9:0, 10:0, 14:0, 16:0, 17:0, 18:1, 18:0 9 -23 8;0 9:0, 14:0, 16:0, 17:0, 18:0 9 -24 8:0, 9:0, 10:0, 14:0, 16:0, 18:0 5, 6, 7, 8, 9, 10 -29 8:0, 9:0, 14:0, 16:0, 18:1, 18:0 6, 7, 8, 9 -30 6:0, 8;0 9:0, 14:0, 15:0, 16:1 16:0, 17:0, 18:1 (tr), 18:0 9, 10 -35 9:0, 12:0, 14:0, 15:0, 16:0, 18:0 7, 9 -41 8;0 9:0, 13:0, 14:0, 15:0, 16:0, 18:1, 18:0, 24:0 6, 7, 8, 9, 10 26:0 (tr), 28:0 (tr) 43 9:0, 10:0, 12:0, 14:0, 15:0, 16:0, 17:0, 18:1, 18:0, 20:0, 24:0 9:0, 14: 0, 15:0, 16:0, 17:0, 18:1, 18:0 9 11 14:0, 15:0, 16:0, 18:1, 18:0 -13 9:0, 14:0, 16:0, 18:1 (tr), 18:0 -17 14:0, 15:0, 16:0, 18:1, 18:0 -20 9:0, 10:0, 12:0, 13:0, 14:2, 14:1, 14:0, 15:0, 16:1, 16:0, 18:1, 18:0, 17:1, 17:0, 18:1, 18:0, 20:1, 20:0, 22:0, 24:0 9 21 9:0, 12:0, 14:0, 15:0, 16:0, 18:1, 18:0 9 (tr) 37 14:0, 16:0, 18:0 -39 9:0, 12:0, 13:0, 14:0, 15:0, 16:0, 18:1, 18:0, 24:0 -68 9:0, 14:0, 16:0, 17:0, 18:1, 18:0 -  Group 1: samples with positive concentration profiles in figure 4B;  Group 2: samples with positive concentration profiles in figure 4D;  Group 3: samples with negative profiles in both of the figures 4B and D.

  of a multidisciplinary research project entitled Human Egyptian Lyon 75 Confluences Mummies (HELYCOM)-Mourir pour renaître, mummification balms of human 76 mummies were analyzed. A total of 47 samples belonging to the Egyptian collection of the 77

129 2 . 4 .

 24 Gas Chromatography coupled to Mass Spectrometry 130 2.4.1. Solvents and reagents 131 All solvents were of the highest purity grade. Hexane, tetrahydrofurane (THF) and N,O-Bis 132 (trimethylsilyl)trifluoroacetamide/Trimethylchlorosilane (BSTFA/TMCS) were supplied by 133 Sigma-Aldrich. Ethanol, dichloromethane (DCM) and diethylether (DEE) were supplied by 134 Merck. 135 136 2.4.2. Samples extraction 137 10 mg of the sample were extracted with 1 mL of dichloromethane (DCM) using ultrasounds 138

  different samples and with different state of conservation (deterioration/damage). 300 301 Acknowledgments 302 The authors want to thank Pr Cathy Vieillescazes and Dr Annie Perraud for their collaboration, 303 the Museum of Confluences in Lyon, the LabEx ARCHIMEDE and Human Egyptian LYon 304 COnfluences Mummies (HELYCOM) -Mourir pour renaître project for the partial funding of 305 this project. A part of this project was supported by LabEx ARCHIMEDE from the 306 Investissement d'Avenir program ANR-11-LABX-0032-01.

Figure 1 :Figure 2 :

 12 Figure 1 : Second extracted spectrum (A) and second associated concentration profile (B) 431 obtained by SIMPLISMA performed on all spectra of mummies samples. 432

Figure 3 :

 3 Figure 3 : Fifth extracted spectrum (A) and fifth associated concentration profile (B) obtained 435 by SIMPLISMA performed on all spectra of mummies samples. 436

Figure 4 :

 4 Figure 4 : Third extracted spectrum (A) and third associated concentration profile (B) and sixth 437 extracted spectrum (C) and sixth associated concentration profile (D) obtained by SIMPLISMA 438 performed on all spectra of mummies samples. 439

Figure 5 :

 5 Figure 5: Principal component analysis (PCA) performed on spectra of 47 mummies samples. 440  Group 1: samples with positive concentration profiles in figure 4B;  Group 2: samples with 441 positive concentration profiles in figure 4D;  Group 3: samples with negative profiles in both 442 of the figures 4B and D. 443

Figure 6 :

 6 Figure 6: Associated loading of PC1 (A) and PC2 (B). 444 445 446

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 6 :

 12346 Figure 1 : Second extracted spectrum (A) and second associated concentration profile (B) 464 obtained by SIMPLISMA performed on all spectra of mummies samples. 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493

Figure 1 :

 1 Figure 1 : Second extracted spectrum (A) and second associated concentration profile (B)obtained by SIMPLISMA performed on all spectra of mummies samples.

Figure 2 :

 2 Figure 2 : Fourth extracted spectrum (A) and fourth associated concentration profile (B) obtained by SIMPLISMA performed on all spectra of mummies samples.

  Figure 2

Figure 3 :

 3 Figure 3 : fifth extracted spectrum (A) and fifth associated concentration profile (B) obtained by SIMPLISMA performed on all spectra of mummies samples.

Figure 4 :

 4 Figure 4 : Third extracted spectrum (A) and third associated concentration profile (B) and sixth extracted spectrum (C) and sixth associated concentration profile (D) obtained by SIMPLISMA performed on all spectra of mummies samples.

  Figure 4

Figure 5 :

 5 Figure 5: Principal component analysis (PCA) performed on spectra of 47 mummies samples.

Figure 6 :

 6 Figure 6: Associated loading of PC1 (A) and PC2 (B).

Table legends 447Table 1 :

 legends1 Presentation of 47 samples of balms from 35 Human mummies. 448

Table 2 :

 2 Fatty matter molecules detected, by gas chromatography and mass spectrometry, on 449 all samples of cluster 1. 450

Table 3 :

 3 Fatty matter molecules detected, by gas chromatography and mass spectrometry, on 451 all samples of cluster 2. 452

Table 1 :

 1 Presentation of 47 samples of balms from 35 Human mummies. 646

	Inventary	N°	Sex	Age	Localisation	Dating 14 C
	30000104 29	Female	Adult	Unknown	Late Period to Ptolemaic Period
	30000107 68	Female	Elderly	Gournah	Unknown
	30000110 14	Male	Adult	Thebes	Graeco-Roman Period
	30000111 34	Male	Elderly	Gournah	Graeco-Roman Period
	30000114	1	Female	Elderly	Thebes	Ptolemaic to Roman Period
	30000125	5	Male	Adult	Gournah	Late Period to Ptolemaic Period
	30000144	6	Male	Mature	Unknown	Ptolemaic to Roman Period
	30000146 35	Female	Adult	Thebes	Unknown
		22
	30000133	23	Male	Adult mature	Gournah	Late Period to Ptolemaic Period
		24
	30000139 15	Female	Young adult	Gournah	Ptolemaic to Roman Period
	30000142 63	Male	Adult	Gournah	Late Period to Ptolemaic Period
	30000143	12 37	Male	Adult	Thebes	Third Intermediate Period to Late Period
	30000148 17	Male	Elderly	Thebes	Unknown
		19
	30000150	20 21	Male	Adult	Deir el-Medineh	Unknown
		33
	30000152	2	Female	Teenager	Thebes	Graeco-Roman Period
	30000286	10 11	Male	Mature	Esneh	Third Intermediate Period to Late Period
	30000306 30	Female	Young adult	Deir el-	Unknown
	30000364 13	Female	Adult	Medineh Unknown	Unknown
	90001402	8 9	Male	Adult	Unknown	Third Intermediate Period to Late Period
	90001635	38 39	Female	Young adult	Unknown	Predynastic
	90001169 41	Male	Young adult	Thebes	Late Period to Ptolemaic Period
	90001259 47	Female	Adult	Kôm Ombo	Ptolemaic to Roman Period
	90001598 53	Female	Adult	Unknown	Nouvel Empire
	90002978A 60	Male	Adult mature	Gournah	Unknown
	90001627 61 Unknown	Child	Rôda	Late Period to Ptolemaic Period
	90001948A	44 45	Male	Adult	Gournah	New Kingdom to Third Intermediate Period
	90001948A 46	Male	Adult	Gournah	New Kingdom to Third Intermediate Period
	90001949A 49	Female	Elderly	Gournah	Graeco-Roman Period
	90001950A 52	Female	Adult	Gournah	Late Period to Ptolemaic Period
	90001951A	50 Unknown 51	Adult	Gournah	Unknown
	90001953A 58	Male	Adult	Gournah	New Kingdom to Third Intermediate Period
	90001956A 56	Female	Young	Gournah	Unknown
	90001957A 57 Unknown	Adult	Gournah	Third Intermediate Period to Late Period
	90001257	43 69	Male	Adult	Gournah	Unknown
	90001597 66	Female	Adult	Unknown	Late Period to Ptolemaic Period
	647
	648
	649

Table 2 :

 2 Fatty matter molecules detected, by gas chromatography and mass spectrometry, 650 on 9 samples of cluster 1. 651

	Ech	Monocarboxylic fatty acid	Dicarboxylic fatty	degraded
			acid	beeswax
	1	8:0, 9:0, 14:0, 15;0, 16:0, 18:1, 18:0	5, 7, 8, 9, 10	-
	2	8:0, 9:0, 10:0, 14:0, 16:0, 18:1, 18:0	4, 5, 6, 7, 8, 9, 10,	-
			11, 12

Table 3 :

 3 Fatty matter molecules detected,by gas chromatography and mass spectrometry, on 656 9 samples of cluster 2. 657

	Ech	Monocarboxylic fatty acid	Dicarboxylic fatty
			acid
	8

Salvado nativitat

Departament d'Enginyeria Qu´ımica EPSEVG, Universitat Polit`ecnica de Catalunya, Av. V´ıctor Balaguer s/n, 08800 Vilanova i la Geltr´u, Barcelona, Spain E-mail address: nativitat.salvado@upc.edu (N. Salvad´o).

List of Three Potential Reviewers

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

8;0 9:0, 10:0, 14:0, 16:0, 17:0, 18:1 (tr), 18:0, 20:0, 22:0, 24:0 5, 6, 7, 8, 9, 10 26:0, 28:0, 30:0 9 8;0 9:0, 10:0, 12:0, 14:0, 15:0, 16:0, 17:0, 18:0 7, 8, 9 24:0, 26:0, 28:0 8:0, 9:0, 10:0, 14:0, 15:0, 16:0, 17:0, 18:1, 18:0 5, 6, 7, 8, 9, 10 -8:0, 9:0, 10:0, 14:0, 15;0, 16:0, 18:1, 18:0, ricinoléique 5, 6, 7, 8, 9, 10 -8:0, 9:0, 10:0, 12:0, 14:0, 15:0, 16:0, 18:0 7, 8, 9, 10 -8;0 9:0, 10:0, 14:0, 16:0, 17:0, 18:1, 18:0 9 -8;0 9:0, 14:0, 16:0, 17:0, 18:0 9 -8:0, 9:0, 10:0, 14:0, 16:0, 18:0 5, 6, 7, 8, 9, 10 -8:0, 9:0, 14:0, 16:0, 18:1, 18:0 6, 7, 8, 9 -6:0, 8;0 9:0, 14:0, 15:0, 16:1 16:0, 17:0, 18:1 (tr), 18:0 9, 10 -9:0, 12:0, 14:0, 15:0, 16:0, 18:0 7, 9 -8;0 9:0, 13:0, 14:0, 15:0, 16:0, 18:1, 18:0, 24:0 6, 7, 8, 9, 10 26:0 (tr), 28:0 (tr) 9:0, 10:0, 12:0, 14:0, 15:0, 16:0, 17:0, 18:1, 18:0, 20:0, 24:0 4, 5, 6, 7, 8, 9, 10 26:0, 28:0 8;0 9:0, 12:0, 13:0, 14:0, 16:0, 18:1, 18:0, 24:0 4, 5, 6, 7, 8, 9 26:0, 28:0, 30:0 9:0, 14:0, 15:0, 16:1, 16:0, 17:0, 18:1, 18:0, 20:1, 20:0, 22:1,

0, 14:0, 15:0, 16:0, 17:0, 18:1 (tr), 18:0 7, 8, 9, 10 -9:0, 12:0, 14:0, 15:0, 16:1, 16:0, 17:0, 18:0, 20:0, 22:0, 24:0 5, 6, 7, 8, 9, 10, 11 -9:0, 14:0, 16:0, 18:1 (tr), 18:0 6, 7, 8, 9, 10 -9:0, 10:0, 12:0, 14:0, 15:0, 16:0, 18:1, 18:0, ricinoléique, 20:0, 22:0, 24:0 5, 6, 7, 8, 9, 10, 11 -9:0, 10:0, 14:0, 15:0, 16:1, 16:0, 17:0, 18:1, 18:0, ricinoléique, 20:0, 22:0, 24:0 5, 6, 7, 8, 9, 10 -9:0, 10:0, 12:0, 14:0, 15:0, 16:0, 18:1 (tr), 18:0, 20:0, 22:0, 24:0 5, 6, 7, 8, 9, 10, 11 -9:0, 12:0, 14:0, 15:0, 16:1, 16:0, 18:1, 18:0, 20:0, 22:0, 24:0 5, 6, 7, 8, 9, 10, 11 -9:0, 14:0, 16:1 (tr), 16:0, 18:0 4, 5, 6, 7, 8, 9, 10, 11 -9:0, 14:0, 15:0, 16:0, 17:0, 18:1, 18:0, 20:0, 22:0, 24:0 5, 6, 7, 8, 9, 10 -

Table 3: Fatty matter molecules detected,by gas chromatography and mass spectrometry, on 9 samples of cluster 2.

Ech Monocarboxylic fatty acid

Dicarboxylic fatty acid 8 9:0, 14: 0, 15:0, 16:0, 17:0, 18:1, 18:0 9 14:0, 15:0, 16:0, 18:1, 18:0 -9:0, 14:0, 16:0, 18:1 (tr), 18:0 -14:0, 15:0, 16:0, 18:1, 18:0 -9:0, 10:0, 12:0, 13:0, 14:2, 14:1, 14:0, 15:0, 16:1, 16:0, 18:1, 18:0, 17:1, 17:0, 18:1, 18:0, 20:1, 20:0, 22:0, 24:0 9 9:0, 12:0, 14:0, 15:0, 16:0, 18:1, 18:0 9 (tr) 14:0, 16:0, 18:0 -9:0, 12:0, 13:0, 14:0, 15:0, 16:0, 18:1, 18:0, 24:0 -9:0, 14:0, 16:0, 17:0, 18:1, 18:0 -

Supplementary Material

Click here to access/download Supplementary Material Supplementary material momies.docx


Fabien Anthelme
email: fabien.anthelme@ird.fr@orcid.org

Ines Carrasquer

Jorge Luis Ceballos

Gwendolyn Peyre

Novel plant communities after glacial retreat in Colombia: (many

Keywords: biogeographic origins, climatic debt, growth forms, Los Nevados national park, post-glacial chronosequence, tropical alpine

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

The environmental constraints applied by arctic and alpine environments on plants at the boundaries of life are severe and multiple [START_REF] Crawford | Plants at the margin: ecological limits and climate change[END_REF]). However, they do not turn necessarily into stress since plants have evolved and developed major adaptations to increase their fitness in this type of ecosystem [START_REF] Larcher | Stress bei pflanzen[END_REF][START_REF] Körner | Alpine plant life: functional plant ecology of high mountain ecosystems 2 nd Ed[END_REF]. With sufficient time available, these species are capable of coping with the effects of natural climatic oscillations through adaptations and migrations, both altitudinal and latitudinal, as proven by the presence of the same taxa in alpine environments of different continents, at different latitudes (e.g., Cerastium; Scheen et al. 2004 and Draba;[START_REF] Jordon-Thaden | Species richness of the globally distributed, arctic-alpine genus Draba L.(Brassicaceae)[END_REF]. However, the Anthropocene era poses a difficult challenge for alpine plants: to adapt to rapidly changing climate and land use or to migrate within an exceptionally short period of time [START_REF] Lenoir | Climate-related range shifts-a global multidimensional synthesis and new research directions[END_REF]. Among the options available to alpine populations at the limits of life in order not to disappear because of global warming, upward migration is the most common scenario (e.g., [START_REF] Steinbauer | Accelerated increase in plant species richness on mountain summits is linked to warming[END_REF][START_REF] Peyre | The fate of páramo plant assemblages in the sky islands of the northern Andes[END_REF]. Indeed, despite having local migration options, e.g., beneficiating of microclimatic refuges proposed by rocks or nurse plants (Anthelme et al. 2014a) the majority of alpine plants are expected to extend their distribution upward, seeking newly available climatic niches, whereas their bottom distribution is negatively affected by thermophilization and competition with species migrating from lower vegetation belts [START_REF] Lenoir | Climate-related range shifts-a global multidimensional synthesis and new research directions[END_REF]. An increasing body of literature evidences a climatic debt, as a result of a significant lag between the pace at which the climate warms and the velocity of plants' upward migration [START_REF] Kuussaari | Extinction debt: a challenge for biodiversity conservation[END_REF][START_REF] Gaüzère | Where do they go? The effects of topography and habitat diversity on reducing climatic debt in birds[END_REF]. For mountain plant species this lag has been split into three types: dispersal lag, establishment lag and extinction lag [START_REF] Alexander | Lags in the response of mountain plant communities to climate change[END_REF]). In the upper alpine zone, where upward migration is assimilated to a primary succession on mostly abiotic terrains [START_REF] Walker | Primary succession and ecosystem rehabilitation[END_REF], this lag has been associated with dispersal limitation (dispersal lag) and a deficit in positive interactions between plants (establishment lag; [START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF][START_REF] Llambí | Vegetation assembly, adaptive strategies and positive interactions during primary succession in the forefield of the last Venezuelan glacier[END_REF]). Nevertheless, the very slow development of alpine soils, especially the slow enrichment in organic matter [START_REF] Khedim | Topsoil organic matter build-up in glacier forelands around the world[END_REF]) means that, even if plants are capable of reaching upward sites, they may not be able to establish because of the lack of organic matter (establishment lag). As a consequence, novel alpine communities at higher elevations are expected to be different from the initial species pool they come from due to (1) the loss of species of the initial pool unable to disperse or establish in novel, deglacierized terrains within a few decades and (2) the intrusion of species that were initially not part of the alpine community but are better adapted to rapid changes (Le [START_REF] Roux | Rapid range expansion and community reorganization in response to warming[END_REF]Peyre et al. 2020). Still, the ecology of páramo plant communities, i.e., the higher tropical alpine belt, close to glaciers, is understudied in comparison with its alpine counterparts outside the tropics [START_REF] Anthelme | Alpine and arctic plant communities: a worldwide perspective[END_REF]. Identifying species' losses and gains is a topical challenge given the high pace of climate warming in this region, which belongs to one of the major hotspots of vascular plant biodiversity in the world: the Tropical Andes [START_REF] Mittermeier | Global biodiversity conservation: the critical role of hotspots[END_REF][START_REF] Peyre | Terrestrial biodiversity hotspots: Challenges and opportunities[END_REF]. From a functional viewpoint, the few studies examining the first decades of plant primary succession after glacial retreat in the alpine tropics evidenced significant changes in species strategies (e.g., from opportunistic to stress tolerant; [START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF][START_REF] Rosero | Multi-taxa colonisation along the foreland of a vanishing equatorial glacier[END_REF] and in growth forms [START_REF] Llambí | Vegetation assembly, adaptive strategies and positive interactions during primary succession in the forefield of the last Venezuelan glacier[END_REF]. This said, the absence of pronounced seasonality and seasonal snowpack as well as the decreasing precipitation at higher elevations in tropical alpine environments have been suggested to result in different functional trait values in the tropics and in temperate-arctic environments [START_REF]Aboveground-trait variations in 11 (sub) alpine plants along a 1000-m elevation gradient in tropical Mexico[END_REF]. The presence of growth forms like large stem/basal rosettes and upright shrubs and the relatively high abundance/diversity of cushion-forming species and large tussock grasses in comparison with similar ecosystems at higher latitudes illustrate these differences [START_REF] Hedberg | Tropical-alpine life-forms of vascular plants[END_REF]Aubert et al. 2014). Some of these tropical alpine growth forms of the Northern Andes are dependent on the presence of highly organic, temporarily water-saturated soils, especially in depressions and relatively flat areas, e.g., some species of the giant stem rosette Espeletia [START_REF] Valencia | Climate vulnerability assessment of the Espeletia complex on Páramo Sky Islands in the Northern Andes[END_REF]) and the cushion-forming Distichia muscoides that forms peatlands [START_REF] Benavides | The influence of climate change on recent peat accumulation patterns of Distichia muscoides cushion bogs in the high-elevation tropical Andes of Colombia[END_REF]. Still, the development of such organic soils is particularly slow, even if possibly accelerated with the effects of warming [START_REF] Benavides | The influence of climate change on recent peat accumulation patterns of Distichia muscoides cushion bogs in the high-elevation tropical Andes of Colombia[END_REF][START_REF] Llambí | Vegetation assembly, adaptive strategies and positive interactions during primary succession in the forefield of the last Venezuelan glacier[END_REF]. Consequently, we expect that the first decades after glacial retreat may not provide adequate soil conditions for the establishment of the tropical alpine growth forms that are the most dependent on highly organic soils.

Another potential driver of primary succession specific to the páramos is their location at a biogeographic crossroad between temperate, tropical and endemic trends. Indeed, in this region, taxa may originate from temperate environments (Holarctic species from the North and Austral-Antarctic species from the South) or from tropical environments, including originating in the tropical páramo itself [START_REF] Sklenář | Tropical and temperate: evolutionary history of páramo flora[END_REF][START_REF] Anthelme | Biogeography of South American highlands[END_REF]. The biogeographic origin of taxa has been observed to influence their elevational distribution in the páramo with, e.g., Holarctic genera becoming more abundant and páramo endemics becoming less abundant at higher elevations [START_REF] Sklenář | Geographic flora elements in the Ecuadorian superpáramo[END_REF]. Taxa of temperate origin may be preadapted to alpine environments, which would make them potentially more efficient for colonizing recently deglacierized terrains than their tropical counterparts, which first diversified in warm, moist environments [START_REF] Sklenář | Tropical and temperate: evolutionary history of páramo flora[END_REF][START_REF] Gehrke | Staying cool: preadaptation to temperate climates required for colonising tropical alpinelike environments[END_REF]Peyre et al. in press; but see [START_REF] Nürk | Are the radiations of temperate lineages in tropical alpine ecosystems pre-adapted?[END_REF]. From this viewpoint, we might expect a deficit in species of tropical biogeographic origin during early primary succession. As a space-for-time substitution approach, the study of plant primary succession along post-glacial chronosequences is a complementary approach to long-term monitoring (e.g., GLORIA project: [START_REF] Pauli | Ecological climate impact research in high mountain environments: GLORIA (global observation research initiative in alpine environments)-its roots, purpose and long-term perspectives[END_REF]. It provides the foundation to assess future trajectories of alpine plant communities and ecosystems [START_REF] Matthews | The ecology of recently deglaciated terrain. A geological approach to glacier forelands and primary succession[END_REF][START_REF] Walker | The use of chronosequences in studies of ecological succession and soil development[END_REF]. This approach has been used frequently and successfully in arctic and temperate alpine regions, on a time step between the Little Ice Age and today, showing a global pattern of increasing biodiversity along with time since glacial retreat (Cauvy-Fraunié and Dangles 2019). The use of chronosequences on shorter time scales, since the acceleration of global warming about 60 years ago, is likely to provide us with valuable information on the capacities of (tropical) alpine plant communities to adapt to abrupt environmental changes at the boundaries of life [START_REF] Cannone | Accelerating climate change impacts on alpine glacier forefield ecosystems in the European Alps[END_REF][START_REF] Anthelme | Living at the edge: increasing stress for plants 2-13 years after the retreat of a tropical glacier[END_REF]. We examined plant primary succession in the most phytodiverse alpine region, worldwide: the north-Andean páramos. In this region, and more generally in the tropical Andes, this approach, which requires close collaboration between glaciologists and ecologists, has been used very little so far, in the Central Andes [START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF], Ecuador [START_REF] Anthelme | Living at the edge: increasing stress for plants 2-13 years after the retreat of a tropical glacier[END_REF][START_REF] Rosero | Multi-taxa colonisation along the foreland of a vanishing equatorial glacier[END_REF]Venezuela (Llambi et al. 2021). Succession was studied along the best-documented post-glacial chronosequence in Colombia, at the foot of glacier Conejeras, located at the slopes of the Santa Isabel volcano [START_REF] Ceballos | Fast shrinkage of tropical glaciers in Colombia[END_REF][START_REF] Rabatel | Toward an imminent extinction of Colombian glaciers? Series A[END_REF]. We focused on six periods: four distributed between 0 and 60 years of deglacierization, within the timeframe of the acceleration of global warming, and two others above and below the frontal moraine of the Little Ice Age (hereafter termed LIA, 169 and nearly 7000 years of deglacierization, respectively; Fig. 1). We aimed at determining to what extent species composition, growth forms and biogeographic origins are distributed along the chronosequence, along with environmental parameters. Our specific objectives were to test to what extent the high velocity of warming, represented by rapid deglacierization, affect positively or negatively the spatial distribution of (1) the plant communities adapted to the vicinity of glaciers, (2) the growth forms characteristics of the alpine tropics, and (3) the proportion of taxa of tropical vs. temperate origins. We tested three specific hypotheses:

(1) Novel plant communities in the recently deglacierized terrain under study are impoverished by the absence of a number of species present in the initial pool, locally (not constrained by time since deglacierization) and enriched by species absent from the initial pool.

(2) The growth forms that require highly organic soils may see their colonization pace slowed down by the establishment lag between the speed of global warming and the slowness of soil formation (climatic debt).

(3) Taxa of temperate biogeographic origin may establish better than taxa of tropical origin recently after deglacierization because they are better pre-adapted to the climate of alpine environments.

MATERIAL AND METHODS

Study area and study site

In Colombia, the northern part of the tropical Andes divides into three main mountain systems: the western, central and eastern cordilleras, above whose treeline the paramo -the local alpine zone-occupies 1% of the national, terrestrial territory [START_REF] Peyre | Mapping the páramo land-cover in the Northern Andes[END_REF]. The Colombian paramo is highly fragmented and consists of a sky island system above 3000-4000 m a.s.l. (Anthelme et al. 2014b;[START_REF] Sklenář | Island biogeography of tropical alpine floras[END_REF][START_REF] Peyre | Mapping the páramo land-cover in the Northern Andes[END_REF].

The vegetation of the Colombian páramo is under the influence of harsh environmental conditions that reduce the diversity of ecological niches for plant species. These include very low minimum temperatures each night, very high daily temperature ranges and heterogeneous, generally nutrient-poor soils [START_REF] Luteyn | Páramos: a checklist of plant diversity, geographical distribution, and botanical literature[END_REF][START_REF] Körner | Alpine plant life: functional plant ecology of high mountain ecosystems 2 nd Ed[END_REF]. The vegetation also faces disturbances related to frequent freeze/thaw cycles and generally steep surfaces. The páramos are also subject to constraints characteristic of tropical alpine environments, such as low annual seasonality, lack of seasonal snow cover, reduced annual precipitation at higher elevations and higher UV radiation [START_REF] Körner | Alpine plant life: functional plant ecology of high mountain ecosystems 2 nd Ed[END_REF][START_REF] Buytaert | Potential impacts of climate change on the environmental services of humid tropical alpine regions[END_REF][START_REF] Vuille | Rapid decline of snow and ice in the tropical Andes-Impacts, uncertainties and challenges ahead[END_REF]. Four glacier complexes, located in the Central and the Eastern cordilleras, remain in Colombia, which are about to disappear within decades because of the effects of accelerating climate warming [START_REF] Rabatel | Toward an imminent extinction of Colombian glaciers? Series A[END_REF]).

The Here, the plant communities belong the midpáramo and the superpáramo vegetation belts (sensu [START_REF] Luteyn | Páramos: a checklist of plant diversity, geographical distribution, and botanical literature[END_REF], which corresponds to the alpine belt in temperate environments. Lower elevations are represented by the upper zonal midpáramo represented by bunch grasslands with the large tussocks Calamagrostis recta (Kunth) Trin. ex Steud. and C. effusa (Kunth) Steud, the upright shrub Diplostephium schultzii Wedd. and the stem rosette Espeletia hartwegiana Cuatrec. At the higher elevations of the superpáramo the vegetation cover is patchy and scarce, represented by various species including the erect herbs Agrostis foliata Hook.f. and Cerastium floccosum Benth and the acaulescent rosette Draba spp. and Erigeron chionophilus Wedd. [START_REF] Salamanca | The paramo vegetation of the volcanic Ruiz-Tolima massif[END_REF]. The sampled site does not include the superpáramo blue grasslands, located between the two above-described communities and characterized by Festuca breviaristata Pilg. and Agrostis araucana Phil. The páramo of Los Nevados harbors the endemic angiosperms Senecio isabelis S. Díaz, Draba pennel-hazenii O.E.Schulz and Lupinus ruizensis C.P.Sm., and the endemic hummingbird Oxypogon stuebelii Meyer.

The Conejeras glacier has undergone two significant mass growth during the Holocene, one 7400-6050 years ago during a climate anomaly and a more recent one during the LIA, around 1850 in Colombia [START_REF] Thouret | Late Quaternary glacial stades in the Cordillera Central, Colombia, based on glacial geomorphology, tephra-soil stratigraphy, palynology, and radiocarbon dating[END_REF][START_REF] Florez | Colombia: evolución de sus relieves y modelados[END_REF]. In 1850, the Santa Isabel ice cap had an area of 27.8 km 2 . It decreased to 10.8 km 2 in 1976. Since the acceleration of global warming in the 1960's, the ice cap has further decreased to only 1 km 2 in 2016, of which less than 0.2 km 2 is represented by the Conejeras glacier. In 2020, the glacier has been split into two fragments (Supplementary Material 2) and it is expected to disappear in the coming years [START_REF] Rabatel | Toward an imminent extinction of Colombian glaciers? Series A[END_REF].

We delimited a precise post-glacial chronosequence focused on the time scale of the acceleration of global warming (since the 60's) and including the LIA moraine. For this purpose, we extracted the glacier limits in 1959,1975,1991,2007 and 2019 using the available aerial photographs at the Agustin Codazzi Geographical Institute. The photographs were georeferenced and the glacier's limited manually created on QGIS (URL: http://qgis.osgeo.org). To compare the recently deglacierized zones with that before and after the LIA (approx. 1850 AD) we used the LIA delimitation, provided by the Institute of Hydrology, Meteorology and Environmental Studies (IDEAM). The resulting chronosequences had six deglacierization zones: 2007-2019 (zone 1), 1991-2007 (zone 2), 1975-1991 (zone 3), 1959-1975 (zone 4), ~1850 (zone 5, just above the LIA moraine) and ~7000 years before present (zone 6, just below the LIA moraine; Fig. 1; Supplementary Material 1). The use of glacial advances to date the onset of primary succession in the different areas was considered relevant because they represent the last major disturbances, with the last locally observed lava flow dating to 5800 years BP and not affecting the glacier-affected area. [START_REF] Kuhry | A paleobotanical and palynological study of Holocene peat from the El Bosque mire, located in a volcanic area of the Cordillera Central of Colombia[END_REF].

Zones 1, 2, 3 and 4 represented recent deglacierization at the time of the acceleration of warming; they were located at similar elevation. Twelve temperature loggers (HOBO Pendant® 8K data logger), two at each zone, were used to described temperature conditions along the chronosequence. They were installed 5 cm below the ground between 22 June 2018 and 12 July 2018 and they recorded temperature every 10 minutes. Data in zones 1-2-3-4 showed that daily minimum and maximum soil temperature increased with time since glacier retreat (e.g., between 0.67°C in zone 1 and 1.24°C in zone 4 for minimum temperature; Supplementary Material 5). Zone 5 and zone 6 represented older zones constrained (5) or not (6) by the presence of a glacier at the time of the LIA. They were located at the same elevation, only separated by the frontal moraine of the little ice age (Supplementary Material 1). Daily minimum temperature was higher in zone 6 (4.31°C vs. 3.75°C in zone 5) whereas daily maximum temperature was higher in zone 5 (10.45°C vs. 8.43°C in zone 6; Supplementary Material 5). The elevational difference between zones 1-2-3-4 and zones 5-6 approximates 350 m, a parameter that might explain part of the variation in our biotic and abiotic measurements along the chronosequence.

Sampling design and data collection

To sample a representative area of each zone, we haphazardly selected 15 squared vegetation plots of 4 m 2 after estimating that this surface was representative of the species richness, locally (see the species-area relationship curve built on plots of 4, 2, 1 and 0.25m 2 in size and showing no strong species enrichment beyond 4m 2 ; in Supplementary Material 3). The minimum distance between plots was 30 m and each plot should contain at least one vascular species, to avoid plots devoid of vegetation in the recently deglacierized zones, which would have weakened statistical analyses. As a result, the species richness in these zones might be slightly overestimated, which is not of significant concern given the objectives and the hypotheses of the study. Five plots in each zone were selected in hydromorphic areas (presence of water as a periodic stream or a pond) as this habitat was frequently observed in the sampling. Areas dominated by flagstones were avoided as being extreme constraints for plant colonization. More generally, we avoided sampling plots with extreme topography, not representative of each zone, to reduce the effects of the topography variable, which was not part of the original hypotheses. The sampling resulted in 90 plots distributed into six zones that includes superpáramo communities (Z1, 2, 3, 4) and the upper bunchgrass mi páramo (Z5 and Z6). At each plot, we made the list of the vascular species by provisionally identifying them on the field and by collecting two specimens of each species for definitive identification at the herbarium of the University Los Andes where all samples were kept (Collecting permit AUR 007-19; see the full list in Supplementary Material 4 and photographs of species in the web/mobile application Plantnet: https://identify.plantnet.org/fr/groups/19892283922). Species names refer to the World Flora Online (http://www.worldfloraonline.org; access 28/11/2021). We visually estimated the relative cover of each species at soil surface, as well as the relative cover of biological soil crust (BSC), rock and bare soil, using a grid of 10*10 cm. We also measured slope degree and exposition with a clinometer (Sunnto Tandem). Plants were grouped into nine growth forms inspired by the classifications by [START_REF] Hedberg | Tropical-alpine life-forms of vascular plants[END_REF] and [START_REF] Ramsay | The growth form composition of plant communities in the Ecuadorian páramos[END_REF]: acaulescent rosette (AR), basal rosette (BR), cushion (CU), erect herb (EH), prostrate herb (PH), prostrate shrub (PS), stem rosette (SR), large tussock (LT) and upright shrub (US; Fig. 2). BR were large rosette-like species with a basal leaf rosette and a large inflorescence that can exceed one meter were grouped into this growth form. SR was represented with only one species: E. hartwegiana. PH and PS had essentially a horizontal structure aboveground. LT took into account the species with a diameter generally exceeding 1m in mature individuals. Species were grouped into CU as soon as they had a dense, hemispheric shape, providing a singular microenvironment. Species were also grouped into biogeographical origins -temperate, tropical and cosmopolite-following the classification of [START_REF] Sklenář | Tropical and temperate: evolutionary history of páramo flora[END_REF] at the genus level, which means each species was grouped into the biogeographical origin of its genus. This was done to test to what extent the biogeographic origin of species conditions their establishment success during primary succession (Supplementary Material 4). To measure changes in the abiotic environments and make correlations with biotic data, soil samples were extracted in ten of the plots selected in each zone, between 0 and 30 cm belowground (or less when unavailable), after removing the organic litter (recognizable pieces of plants) when present. Each sample was a combination of five subsamples collected at the corners and the center of each plot (total: approximately 500g). Soils were brought to the Environmental Engineering laboratory of the Los Andes University (Bogota, Colombia) for measurements of pH and soil organic matter content (SOM; %) following the protocol of [START_REF] Jaramillo | Introducción a la ciencia del suelo[END_REF]. In addition, at each plot we measured the soil water content (SWC; %) with vertical probes reaching 7.5 cm belowground (Fieldscout TDR-100; SWC at each plot was the average value of five measurements made at the corners and at the center). Soil and SWC collections were made during two days of stable climatic conditions without significant precipitations. Soil and SWC were measured within this time period. Averaging the two values in each zone showed that minimum temperature increased monotonously between zone 1 (0.67°C) and zone 6 (4.31°C; Supplementary Material 5). Mean temperature also increased between zone 1 (1.78) and zone 6 (6.21) but reached a maximum value in zone 5 (6.70°C).

Data analysis

To test hypothesis 1, changes in plant communities along the chronosequence were assessed with a partition around medoids (PAM) clustering analysis, which establishes clusters, or groups, of sampling units (i.e., vegetation plots) with similar variable values (i.e., species) and distribution (i.e., relative cover) [START_REF] Kaufman | Finding groups in data: an introduction to cluster analysis[END_REF]. We carried out the PAM analysis on a Bray Curtis distance matrix of our plot data with partitions increasing cluster numbers between 1 and 60 and 100 iterations each (package "cluster"). We used the within sum of squares (wss) metric to establish the best partition to divide the dataset. Each resulting cluster, revolving around a plot-defined medoid, was therefore assimilated as a plant community that we described by its characteristic flora. To do so, we calculated the Indval index to establish the probability for a species to be preferentially found in a cluster in comparison to other clusters (R package indicspecies). To consider a species characteristic of a specific group, we put a 0.5 threshold for the IndVal value and ensured that it was substantially less represented in other groups (min 30% less). We refrained from using the term "diagnostic" to qualify the characteristic species because of the local scale of our analysis, however, we discussed our results relating them to broader scale studies on the páramo vegetation and its characteristic species (e.g., [START_REF] Salamanca | The paramo vegetation of the volcanic Ruiz-Tolima massif[END_REF][START_REF] Peyre | Phytoregionalisation of the Andean páramo[END_REF]). To test hypothesis 2, changes in growth forms along the chronosequence and their correlation with environmental variables were evaluated with a principal component analysis on a matrix crossing plots with the relative cover of growth forms (eight variables) and integrating the variables pH, SOM, SWC and BSC (four environmental variables; R package Factoextra). Lacking values of SOM and SWC (five plots out of 15 in each zone) were filled with average values at each zone. Additionally, linear regressions were used to test the correlations between the cover of each growth and the relative content of soil organic matter. To test hypothesis 3, the relationship between the biogeographic origin of taxa and time since deglacierization was evaluated by comparing the relative species richness of temperate, tropical and cosmopolite origins, in each plot and compared between zones. First, we calculated origin frequency within the plot data and rescaled these frequencies on a 0-1 interval. Then, we conducted a Dirichlet regression on the rescaled frequency data in function of the "zone" variable, using the mean age value of each zone: (6,20,36,52,169 and 7000 years;Supplementary Material 6). We used the alternative parametrization, which is often privileged for ecological data (R package DirichletReg; [START_REF] Maier | DirichletReg: Dirichlet Regression for Compositional Data in R[END_REF][START_REF] Douma | Analysing continuous proportions in ecology and evolution: A practical introduction to beta and Dirichlet regression[END_REF]. Finally, we used the Akaike information criterion (AIC) to indicate model performance and predicted frequency values for each origin category along the chronosequence. The data and predicted values were plotted altogether to visually compare changes in plant community composition in term of temperate, tropical and cosmopolitan frequencies along the chronosequence.

RESULTS

Species composition along the chronosequence

A total of 68 vascular species were recorded in the plots along the chronosequence, representing eight growth forms (Fig. 2; Supplementary Material 4). 46 species were of temperate origin, 21 species of tropical origin and one species of cosmopolitan origin (Eryngium humile). Espeletia hartwegiana, Calamagrostis effusa and Senecio isabelis were present only in Z6, below the LIA, among others. The non-native Rumex acetosella (Z1, 5, 6), Bromus lanatus and Taraxacum officinale (Z6) were also present in the sampling. Based on the total within sum of squares, the PAM analysis identified four cluster partitions of species along the chronosequence (Fig. 3). Group 1 comprised 100% of the plots in Z1 and was present in Z2, Z3 and Z4. It had no characteristic species. Group 2 was present in Z2, Z3 and included the majority of the plots of Z4 ( 12). It was characterized by four early-colonizing species: the erect herbs Cerastium floccosum and Agrostis foliata (IndVal: 0.80 and 0.78, respectively) and the acaulescent rosette Erigeron chionophilus (Asteraceae; 0.77). Group 3 gathered all the plots of zone 5 and two plots in Z6. It was characterized by 13 species, notably the prostrate shrubs Pernettya prostrata and Baccharis caespitosa (0.84 and 0.74, respectively) the erect herbs Agrostis haenkeana and Cerastium kunthii, as well as the prostrate herbs Lupinus microphyllus and Lachemilla nivalis (0.84, 0.54, 0.76 and 0.70). Group 4 was present only in Z6. It was characterized by 11 species, among which the large tussock grasses Calamagrostis effusa and Festuca asplundii (1.00 and 0.54), the upright shrubs Diplostephium Schultzii and Pentacalia vernicosa (0.62 and 0.54) and the endemic basal rosette Senecio isabelis (0.62).

Growth forms

Regarding soil parameters, the PCA's axis 1 represented a gradient of time since glacial retreat from the right to the left where mineral, inorganic soils with high abundance of biological soil crust in recently deglacierized terrains were gradually replaced by organic soils with higher water content and a more acidic pH (Fig. 4). Axis 2 separated mainly plots of Z5, in the lower part of axis 2, and plots of Z6, in the upper part of axis 2 and associated with highly organic soils. Along axis 1, all growth forms were associated with the older deglacierization zones on the left of the graph. The growth forms erect herb, prostrate herb, prostrate shrub and, to a lesser extent, cushion were located preferentially in Z5, 169 years after deglacierization. In contrast, the growth forms large tussock, upright shrub and, to a lesser extent, stem rosette were preferentially observed in Z6, in areas deglacierized approximately ~7000 years ago. The relative cover of the two most characteristic growth forms of Z6, large tussock and upright shrub, was positively correlated, linearly, with SOM (R 2 = 0.61 and 0.42, respectively; Fig. 5).

Fig. 4 PCA biplot combining plot (symbols, N = 90), grouped by zone (color) with centroids represented as larger symbols, with growth forms (black, regular) and environmental variables (rock, pH, biological soil crust -BSC-, soil water content -SWC-and soil organic matter -SOM; gray, italic). Variables in black with decreasing length according to their contribution. Full names of growth forms available in Fig. 2 Fig. 5 Relations between the percentage of soil organic matter and the relative cover of the two most characteristic growth forms of zone 6: large tussock and upright shrub. (Red lines: significant linear regressions. SR: regression not applicable because too few values different from 0

Biogeographic origin

The Dirichlet model (Precision estimate: 1.5668; p-value: <2e-16; AIC: -731.7; Supplementary Material 6) predicted a relative increase of the temperate element and a relative decrease of the tropical element along the chronosequence (Fig. 6). In fact, 91% of the recorded species in Z1 plots were of temperate origin on average, while the remaining 8% were of tropical origin (Supplementary Material 4). In Z6, temperate species accounted for 72% of the sampled flora on average, while tropical species encompassed 22% and cosmopolitan species 6%. Temperate species were generally diverse in all zones, across plant families and growth forms. In contrast, tropical species were overall less represented and reached the highest elevations by the means of few species only, such as Belloa kunthiana and Xenophyllum humile (Asteraceae). The frequency trends observed were, in addition, variable among the three categories: tropical species decreased more rapidly than temperate species increased. Finally, cosmopolitan species were little represented, and followed a more linear decreasing trend along the chronosequence.

Fig. 6

Compositional changes at the biogeographic origin level along the chronosequence. Progressive changes in frequency of temperate, tropical and cosmopolitan species from the vegetation plot data and Dirichlet regression results with age of deglacierization

DISCUSSION

Collaborations between glaciologists and ecologists allow the use of a powerful space-for-time substitution approach: post-glacial chronosequences [START_REF] Matthews | The ecology of recently deglaciated terrain. A geological approach to glacier forelands and primary succession[END_REF]. Applied here to the highly phytodiverse páramo, it provides valuable insights on the climatic debt endured by plant communities during primary succession after rapid deglacierization, which will itself affect provisional, regulation and cultural ecosystem services in downstream areas [START_REF] Luckman | Darkening peaks: glacier retreat, science, and society[END_REF]. We used this method to test to what extent the novel tropical alpine communities establishing after the retreat of glaciers differ from initial communities in terms of species, growth forms and biogeographic origins. However, it is to be noted that, despite our sampling efforts to reduce their effects, the parameters local topography, elevation and geomorphology can still explain a significant part of the observed variations.

Novel plant communities after deglacierization

Comparing novel plant communities that established less than 169 years ago with the (super)páramo communities described locally in the literature [START_REF] Salamanca | The paramo vegetation of the volcanic Ruiz-Tolima massif[END_REF][START_REF] Peyre | Phytoregionalisation of the Andean páramo[END_REF]) is a relevant way to determine how time constrains the establishment of tropical alpine communities after deglacierization. The fact that, above the LIA moraine, the first colonizing species Agrostis foliata, Cerastium floccosum, Culcitium canescens, Erigeron chionophilus, (group 2), Baccharis caespitosa, Hypochaeris sessiliflora, Lachemilla nivalis, Luzula racemosa (group 3), were characteristic species both in our sampling and in the above-cited works shows, to some extent, that the novel communities establishing after deglacierization develop common features with pre-exisiting communities. In contrast, the absence above the LIA moraine of Calamagrostis effusa (characteristic of the upper bunchgrass páramo; [START_REF] Salamanca | The paramo vegetation of the volcanic Ruiz-Tolima massif[END_REF], Draba pennell-hazenii (superpáramo endemic of the Central Cordillera) and Senecio isabelis (endemic of the páramo of los Nevados) reveals significant losses in the novel communities established less than 169 years ago in comparison with the initial species pool. In parallel, the presence of Rumex acetosella, a non-native species from Europe, in the most recently deglacierized zone of the sampling is a flagship example of the rapid establishment of a few, novel species in the tropical alpine communities. R. acetosella has been shown to be one of the most efficient invaders of tropical alpine communities in the tropical high Andes, with negative, competitive effects on native species [START_REF] Llambí | Reciprocal interactions between a facilitator, natives, and exotics in tropical alpine plant communities[END_REF]2020). Its presence shortly after deglacierization underlines its high capacity to rapidly colonize disturbed areas and contribute generating non-analogous tropical alpine communities. As well, the presence of the non-native Bromus lanatus in Z5 confirm the capacity for exotic plants to colonize alpine environments (see also [START_REF] Cavieres | The role of plant-plant facilitation in non-native plant invasions[END_REF]. It should be noted, however, that R. acetosella was not observed along recent chronosequences in Ecuador and Venezuela [START_REF] Anthelme | Living at the edge: increasing stress for plants 2-13 years after the retreat of a tropical glacier[END_REF][START_REF] Llambí | Vegetation assembly, adaptive strategies and positive interactions during primary succession in the forefield of the last Venezuelan glacier[END_REF], although being frequently observed in these regions [START_REF] Cuesta | Latitudinal and altitudinal patterns of plant community diversity on mountain summits across the tropical Andes[END_REF]. These latter chronosequences are far less visited by tourists than the Colombian chronosequence under study (pers. obs.). The relationship tourisminvasive plants may thus deserve a focus in further investigations of this type. Accordingly, from a taxonomical viewpoint and partially in line with hypothesis 1, the novel tropical alpine communities established in less than 169 years are not the same as the initial plant communities observed in the same superpáramo but not constrained by time (see also [START_REF] Losapio | The consequences of glacier retreat are uneven between plant species[END_REF]. Whereas species losses include characteristic, endemic and iconic species, whose absence may be detrimental to the completeness of the community, the few species gains concern only non-native species, with potential detrimental effects on the community organization. Further research on the species that are not able to colonize the chronosequence are needed, especially if they are present or absent at similar elevation but outside the chronosequence, in order to evaluate more precisely their future distribution in relation with their incapacity to rapidly colonize deglacierized areas. From a methodological viewpoint, it is to be noted that the elevational difference between zone 4 and zone 5 might explain in part the variations found, underlining the potential limitation of using this specific chronosequence. However, the gradual changes observed in soil temperature plant community assembly along the chronosequence (see Supplementary Material 5 and Fig. 4) suggests that this parameter is not of major influence on our data.

Establishment lag for tropical alpine growth forms

Primary succession after rapid deglacierization generates differential colonization velocities among growth forms, primarily because of a dispersal lag, giving priority to wind-dispersed species [START_REF] Alexander | Lags in the response of mountain plant communities to climate change[END_REF][START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF]. The fact that the growth forms upright shrub and large tussock, known as characteristic of tropical alpine regions [START_REF] Hedberg | Tropical-alpine life-forms of vascular plants[END_REF][START_REF] Rada | Plant functional diversity in tropical Andean páramos[END_REF], were found preferentially before the LIA and were observed in plots with higher soil organic matter content indicates, apart from dispersal constraints, that tropical alpine growth forms may be less able than more generalist growth forms to perform a rapid succession after deglacierization. Despite insufficient presence of the growth form stem rosette in our sampled plots to provide a significant statistical test, we expect that this tropical alpine growth form might as well require high levels of soil organic matter to develop (see also [START_REF] Urbina | Simulated small scale disturbances increase decomposition rates and facilitates invasive species encroachment in a high elevation tropical Andean peatland[END_REF]. Organic soil formation in the tropical high Andes during primary succession is a very slow process, which requires thousands of years [START_REF] Benavides | The influence of climate change on recent peat accumulation patterns of Distichia muscoides cushion bogs in the high-elevation tropical Andes of Colombia[END_REF][START_REF] Khedim | Topsoil organic matter build-up in glacier forelands around the world[END_REF]. Our data suggest that the soils above the LIA moraine (169 years and less) were not sufficiently developed to allow the frequent presence of dependent growth forms to establish themselves despite having the capacity to disperse easily from individuals located on the other side of the moraine (see also [START_REF] Llambí | Vegetation assembly, adaptive strategies and positive interactions during primary succession in the forefield of the last Venezuelan glacier[END_REF]. Accordingly, in line with hypothesis 2, there would be a long-term constraint for plant establishment (establishment lag; [START_REF] Alexander | Lags in the response of mountain plant communities to climate change[END_REF]) and this constraint would be growth form-specific, in accordance with results observed in nearby Ecuador [START_REF] Sklenář | Distribution changes in páramo plants from the equatorial high Andes in response to increasing temperature and humidity variation since 1880[END_REF]. It may be exacerbated by the fact that these large, tropical alpine growth forms have been shown to increase the content in soil organic matter, and associated reduction in pH (see also [START_REF] Llambí | Vegetation assembly, adaptive strategies and positive interactions during primary succession in the forefield of the last Venezuelan glacier[END_REF], thus facilitating potentially the establishment of species requiring this type of habitat [START_REF] Mora | Giant stem rosettes have strong facilitation effects on alpine plant communities in the tropical Andes[END_REF][START_REF] Llambí | Interactions between nurse-plants and an exotic invader along a tropical alpine elevation gradient: growth-form matters[END_REF]. Further soil analyses along chronosequences are expected to add consistence to this exploratory result. It should be remembered, however, that in addition to the effect of soil on plant establishment, other environmental parameters influence the success of colonization of new alpine plant communities in the tropics. These include variations in dispersal abilities, cold and drought tolerance/avoidance, and the presence or absence of nurse plants to mitigate the effects of extreme local climatic conditions [START_REF] Squeo | Freezing tolerance and avoidance in high tropical Andean plants: is it equally represented in species with different plant height?[END_REF][START_REF] Zimmer | Time lag between glacial retreat and upward migration alters tropical alpine communities[END_REF][START_REF] Rada | Plant functional diversity in tropical Andean páramos[END_REF].

Biogeography influences the establishment success

Taxa of tropical origin were under-represented compared to temperate species at the beginning of the succession and this imbalance was progressively corrected along the chronosequence. It shows that the first group may face a greater climatic debt than temperate species along the Conejeras chronosequence. This interpretation is consistent with biogeographic analyses showing that temperate alpine taxa, genetically preadapted to alpine climates, are better candidates for the colonization of tropical alpine environments after deglacierization [START_REF] Sklenář | Tropical and temperate: evolutionary history of páramo flora[END_REF][START_REF] Gehrke | Staying cool: preadaptation to temperate climates required for colonising tropical alpinelike environments[END_REF][START_REF]Aboveground-trait variations in 11 (sub) alpine plants along a 1000-m elevation gradient in tropical Mexico[END_REF]. Our results moreover support previous observations that the temperate vs. tropical element ratio increases from the midpáramo to the superpáramo belt, recorded from 50/50 to 70/30 for Ecuador [START_REF] Sklenář | Tropical and temperate: evolutionary history of páramo flora[END_REF], even though our numbers suggest a more drastic drop for tropical taxa in Colombia. Additional records will help generalize the pattern and identify if historical factors, such as past corridors with Central America, help explaining the ratio difference along the latitudinal gradient. In fact, a recent study on the climatic refugia of páramo plant species at the regional scale has shown that tropical taxa likely found more refugia in the southern and central páramos of Peru and Ecuador, while taxa or temperate origin showed more diverse refugia patterns throughout the region, with an important refugium in the vicinity of the Nevados National Park (Peyre et al. in press). Finally, the absence of Espeletia hartwegiana and Diplostephium schultzii above the limit of the LIA may explain (1) because of their tropical origin, which would make them less tolerant to elevation stress (see absence of Espeletia spp. in the highest peaks of the long-term observatories GLORIA in the northern Andes; [START_REF] Cuesta | Latitudinal and altitudinal patterns of plant community diversity on mountain summits across the tropical Andes[END_REF]) and (2) because they develop large, tropical alpine growth forms that require highly organic soils (see section 4.2. and Fig. 2). For these reasons, E. hartwegiana and D. schultzii, two species iconic of the paramo of Los Nevados, can be considered as sentinel species when studying the effects of warming on the future distribution of tropical alpine plant communities.

Conclusion

According to the Colombian chronosequence, the primary succession initiated by the rapid retreat of the Conejeras glacier, consecutive to rapid warming, generates novel plant communities that differ from initial tropical alpine communities in terms of species and growth form composition and in terms of biogeographic origins. Especially, the establishment lag generated by the slow enrichment of soil organic matter seems to affect negatively and preferentially (1) characteristic, endemic and iconic species, (2) the tropical growth forms upright shrub and large tussock, and (3) taxa of tropical biogeographic origin. Therefore, and if we assume that the chronosequence used is a pertinent space-for-time substitution, the novel plant communities, which represent the future of the tropical alpine communities, could face substantial taxonomic and functional losses because of an exacerbated climatic debt. In return, they receive only a few, questionable quantitative gains, through the establishment of non-native species in the novel plant communities. This is illustrated by a picture taken in November 2021 where the first vascular species to colonize the glacier forefront, one year after glacial retreat, is Poa annua L., a non-native species (pers. obs.; SI2).

Fig. 1

 1 Fig. 1 Post-glacial chronosequence of the rapidly retreating Conejeras glacier, divided into six deglacierization zones, from 2019 to ~7000 years BP (below the LIA moraine in red: zone 6)

  selected study area is the National Park of Los Nevados, located in the Cordillera Central of Colombia and which harbors three of the few remaining glaciers in Colombia (Rabatel et al. 2018). Our study site is the páramo located on the western side of the Volcan Santa Isabel (4965 m) between 4312 and 4716 m a.s.l. (Supplementary Material 1), within the area deglacierized between the LIA and current time below the Conejeras Glacier, approximately 25 km south-east of the city of Manizales (4° 49' 12.9324'' N 75° 22' 33.2112'' W). Mean annual temperature at the foot of the glacier reaches 1.3°C with very little seasonal variation (see also our data in the next section) whereas annual precipitations reach 1025 mm. During the dry season (December-January and June-August) monthly precipitations reach 75 mm per month, whereas it can exceed 150 mm during the rest of the year (wet season; Morán-Tejeda et al. 2018).

Fig. 2

 2 Fig. 2 Nine growth forms observed along the post-glacial chronosequence, with examples. A: acaulescent rosette (Hypochaeris sessiliflora); B: basal rosette (Lupinus alopecuroides); C: cushion (Plantago rigida); D: erect herb (Senecio latiflorus); E: large tussock (Calamagrostis effusa); F: prostrate herb (Geranium sibbaldioides); G: prostrate shrub (Pernettya prostrata); H: stem rosette (Espeletia hartwegiana); I: upright shrub (Diplostephium schultzii). Photographs of all species available at: https://identify.plantnet.org/fr/groups/19892283922

Fig. 3

 3 Fig. 3 Distribution of four plant communities in each of the deglacierization zones, as a result of the PAM clustering analysis

This study does not involve research on human participants or animals.Informed consentInformed consent was obtained from all individual participants included in the study

Acknowledgments

This work was funded by the Joint International laboratory LMI BIOINCA, the University of Los Andes and the IRD. The authors thank the National park Los Nevados for providing access to the study sites.

Statements and declaration

Data availability statement All relevant data may be accessible via the VegAndes database upon request. Photographs of the species taken into consideration in the analysis are available at https://identify.plantnet.org/fr/groups/19892283922

Conflict of Interests/Competing Interests

The authors declare that they have no conflict of interest/competing interests.

Ethical approval

Additional supporting information in the online version of this article (see "Supplementary Material") contains the following: ESM_1: description of the deglacierization zones, ESM_2 -Photograph of the glacier forefront, ESM_3: Species-area relationship, ESM_4: Species list and species characteristics, ESM_5: Temperature at each zone. ESM_6: Details of the Dirichlet regression presented in Fig. 6