Jérémy Migliore
email: jeremy.migliore@imbe.fr

Alex Baumel

Agathe Leriche

Marianick Juin

Frederic Medail

Frédéric Médail

Surviving glaciations in the Mediterranean region: an alternative to the long-term refugia hypothesis

Keywords:

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Mediterranean-climate ecosystems harbour the world's richest extra-tropical floras, but contemporary ecological factors and environmental heterogeneity do not entirely explain their remarkable biodiversity and endemism or their spatial patterns [START_REF] Rundel | Mediterranean biomes: evolution of their vegetation, floras, and climate[END_REF]. It has been suggested for the Cape Region and south-western Australia that long-term environmental stability allowed ancient lineages to continue their diversification [START_REF] Cowling | Variation in plant diversity in mediterraneanclimate ecosystems: the role of climatic and topographical stability[END_REF]. By contrast, clades are generally younger in the Mediterranean region, leading to the hypothesis that environmental instability has prevented the long-term accumulation of diversity [START_REF] Valente | Contrasting evolutionary hypotheses between two mediterraneanclimate floristic hotspots: the Cape of southern Africa and the Mediterranean Basin[END_REF][START_REF] Cowling | Variation in plant diversity in mediterraneanclimate ecosystems: the role of climatic and topographical stability[END_REF]. Indeed, several old lineages, i.e. those of early to middle Cenozoic origins, were pruned by extinction around 3.2 Ma, during the transition from a tropical to a Mediterranean climate [START_REF] Mijarra | Floristic changes in the Iberian Peninsula and Balearic Islands (south-west Europe) during the Cenozoic[END_REF][START_REF] Rodríguez-Sánchez | Past tree range dynamics in the Iberian Peninsula inferred through phylogeography and palaeodistribution modelling: A review[END_REF]. The Quaternary glaciations starting at ca. 2.7 Ma subsequently led to an increased climatic severity characterized by more arid conditions and a mean temperature of the coldest month about 12°C cooler than at present [START_REF] Wu | Climatic changes in Eurasia and Africa at the last glacial maximum and mid-Holocene: reconstruction from pollen data using inverse vegetation modelling[END_REF]. These drastic climatic conditions have favoured open landscapes of grasslands and steppe vegetation to the detriment of more thermophilous and woody vegetation [START_REF] Quézel | Ecologie et biogéographie des forêts du bassin méditerranéen[END_REF][START_REF] Thompson | Plant evolution in the Mediterranean[END_REF]. Nevertheless, recent phylogenetic evidence suggest that most of plant clades of the Mediterranean region diversified between the Miocene and Pliocene (from 23 Ma), i.e. before the onset of the Mediterranean climate [START_REF] Vargas | Phylogenetic evidence for a Miocene origin of Mediterranean lineages: species diversity, reproductive traits and geographical isolation[END_REF], so that extinction of older lineages was counteracted by diversifying processes in younger lineages.

Species Distribution Models (SDM) fitted with palaeoclimatic explanatory variables can provide insights into the processes shaping modern genetic and species diversity while tracking the contraction/expansion of populations in response to successive Quaternary glacial/interglacial oscillations [START_REF] Waltari | Locating Pleistocene refugia: comparing phylogeographic and ecological niche model predictions[END_REF][START_REF] Carnaval | Stability predicts genetic diversity in the Brazilian Atlantic Forest Hotspot[END_REF]. The hypothesis of long-term refugia proposes the stability of suitable areas that allow viable populations to persist through climatic oscillations [START_REF] Ashcroft | Identifying refugia from climate change[END_REF][START_REF] Stewart | Refugia revisited: individualistic responses of species in space and time[END_REF]. Within the framework of refuge theory, we should expect to find a positive correlation between the stability of suitable areas and genetic singularity, since long term populations tend to harbour more diversity than recently expanding ones [START_REF] Svenning | Applications of species distribution modeling to paleobiology[END_REF][START_REF] Hampe | Climate refugia: from the Last Glacial Maximum to the twenty-first century[END_REF][START_REF] Gavin | Climate refugia: joint inference from fossil records, species distribution models and phylogeography[END_REF]. These long-term refugia have been defined as 'phylogeographical hotspots', i.e. significant reservoirs of unique genetic diversity of Mediterranean plant species that probably played a key role in speciation processes [START_REF] Médail | Glacial refugia influence plant diversity patterns in the Mediterranean Basin[END_REF].

The complex spatio-temporal trajectories of species and populations in and out of refugia need to be better integrated [START_REF] Gavin | Climate refugia: joint inference from fossil records, species distribution models and phylogeography[END_REF]. Combining palaeoclimatic SDM and phylogeography facilitates inferring the existence of putative Quaternary long-term refugia. This was the case for the olive tree (Olea europaea L.), of which a large part of its populations located in the south-western and the coastal eastern Mediterranean persisted under continuously suitable areas during the Last Inter-Glacial (LIG, ca. 116-130 ka BP), the Last Glacial Maximum (LGM, ca. 19-26.5 ka), and current climatic conditions (Besnard et al., 2013).

The simultaneous analysis of both SDM and genetic imprints left by range dynamics have been rarely performed on the scale of the entire Mediterranean and even less concerning widespread thermophilous plants. In this context, we have chosen to focus on the myrtle, or Myrtus communis L. (Myrtaceae), because it is one of the few Mediterranean thermophilous evergreen shrubs occurring in coastal and lowland environments, below an altitude of 500 m. Indeed, a few common woody plants of the thermo-Mediterranean shrublands have been studied using a phylogeographical approach (Pinus pinea L.: Vendramin et al., 2008;Cistus monspeliensis L.: Fernandez-Mazuecos & Vargas, 2010; Chamaerops humilis L.: García-Castaño et al., 2014;[START_REF] Guzmán | The colonization history of the Mediterranean dwarf palm (Chamaerops humilis L., Palmae)[END_REF]. The myrtle is characterized by a long biogeographical history, having persisted in the Mediterranean region through the Cenozoic since at least the Neogene period (Migliore et al., 2012). Found primarily in low-altitude habitats, this shrub has probably been highly impacted both by climatic and habitat changes through latitudinal shifts of climatic conditions and eustatic oscillations of the Mediterranean sea [START_REF] Pirazzoli | A review of possible eustatic, isostatic and tectonic contributions in eight late-Holocene relative sea-level histories from the Mediterranean area[END_REF]. Thus, the framework of the expansion/contraction model [START_REF] Bennett | What do we mean by 'refugia'?[END_REF][START_REF] Gavin | Climate refugia: joint inference from fossil records, species distribution models and phylogeography[END_REF] could be relevant for past Myrtus range dynamics, especially during the Quaternary climatic oscillations.

To examine the biogeographical responses of this Mediterranean shrub to severe Quaternary climate conditions, we used MAXENT species distribution models, and detected the current and past suitable areas for M. communis on the scale of its entire Mediterranean and Macaronesian distribution range. From this we inferred potential range shifts and identified refugia on the basis of bioclimatic suitability during the Last Inter-Glacial, the Last Glacial Maximum, and the Mid-Holocene (ca. 6 ka BP). We compared the implications of the SDM study with the phylogeography of M. communis based on plastid and nuclear DNA sequences as well as multilocus AFLP (Amplified Fragment Length Polymorphism) markers.

MATERIAL AND METHODS

MODELLING OF SUITABLE AREAS FOR MYRTUS COMMUNIS

Occurrence data

The study area was divided into a grid of 3,019 1 km × 1 km cells using ARCGIS 10.2 (ESRI, Redlands, USA). All the cells under consideration were at least 75% emerged land area. An occurrence dataset, consisting of 4,168 observations of M. communis, was (Figure 1A), collected from the following sources: the Global Biodiversity Information Facility data portal (2,011 occurrences; http://www.gbif.org/), and our own field GPS data (179 Mediterranean occurrences), GPS data from the French CBNMED database (758 occurrences; http://flore.silene.eu), the Atlas of the Aegean Flora (462 occurrences; A. Strid, unpubl. data), Flora Croatica (660 occurrences; [START_REF] Nikolić | Myrtus communis distribution in Croatia[END_REF], and scientific literature on the flora of northern Africa and Italy (98 occurrences). For modelling purposes, the dataset was reduced to one occurrence per unit cell (n = 3019).

Bioclimatic variables

Six variables related to the ecological requirements of the species were chosen and used to fit the bioclimatic distribution model of M. communis. Isothermality (BIO3) reflects the influence of larger or smaller temperature fluctuations within a month, relative to the year (generally useful for insular and maritime environments). Minimal temperature of the coldest month (BIO6) quantifies potentially lethal frost events and generalized stress due to low temperatures. Annual temperature range (BIO7) is useful for examining the effect of the ranges of annual temperature extremes. Annual precipitation (BIO12) approximates average water availability, and precipitation of the driest month (BIO14) describes the extremes associated with potentially lethal drought events and quantifies stress due to low water availability. Finally, because M. communis requires warmer temperatures during the wettest period of year for its germination, mean temperature in the wettest quarter (BIO8) was also considered. These bioclimatic variables were extracted from the Bioclim dataset, provided by WORLDCLIM 1.4 in a GIS-based raster format (1 x 1 km resolution; [START_REF] Hijmans | Very High resolution interpolated climate surfaces for global land areas[END_REF][START_REF] O'donnell | USGS Data Series 691: Bioclimatic Predictors for Supporting Ecological Applications in the Conterminous United States[END_REF]. As a safeguard against redundancy among our variables, we checked that all pairwise Pearson correlation coefficients between any two of the six selected variables were below 0.71.

Modelling method

Maximum Entropy algorithm as implemented in MAXENT 3.3.1 was used to define suitable areas for M. communis [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF][START_REF] Phillips | Modeling of species distributions with MAXENT: new extensions and a comprehensive evaluation[END_REF]. MAXENT estimates species' suitable areas by finding the distribution of maximum entropy (i.e. closest to uniform) subject to the constraint that the expected value of each environmental variable (or its transform and/or interactions) under the estimated distribution matches its empirical average. The method has been shown to be well adapted to presence-only data and has consistently demonstrated performance that is competitive with other methods [START_REF] Elith | Novel methods improve prediction of species' distributions from occurrence data[END_REF]but see Yackulic et al., 2013;[START_REF] Phillips | Opening the black box: an opensource release of MAXENT[END_REF]. Several models were fitted by investigating all feature types (linear, quadratic, and product) and all combinations of features, except for the threshold feature, which is known to fit overly conservative models. Ten replicates were performed for each model (crossvalidation method). The machine learning fitting process was set to end when training gain fell below the threshold value of 0.0001, or after 3,000 iterations. The user-defined number of background points was set to 50,000, and individual points were chosen to correct for latitudinal and sampling biases (i.e. an oriented random selection of 50,000 points across the entire study area, paying attention to cell area, which is a function of latitude, and to sampling bias). The extrapolate function was set to 'no'.

Model selection

The Akaike Information Criterion (AIC) was calculated for all models and averaged by combination to define the most parsimonious combination of features. Area Under the Curve (AUC) and standard deviation metrics were checked for congruence with AIC results. The final model of the current suitable areas for M. communis was fitted using the selected combination of features, with 50 replicates (identical parameters to the fitting procedure). A Jackknife analysis was performed to evaluate variable contributions to the model, and response curves were created to assess how each variable influences the level of suitability.

Model projections and persistence assessment

The final model was projected using prior climate history to determine the suitable areas for M. communis during the Mid-Holocene, LGM and LIG periods. It was assumed that the ecological requirements of M. communis have remained similar over the last Pleistocene climatic cycles [START_REF] Valiente-Banuet | Modern Quaternary plant lineages promote diversity through facilitation of ancient Tertiary lineages[END_REF][START_REF] Nogués-Bravo | Predicting the past distribution of species climatic niches[END_REF]. Reconstructed bioclimatic variables were provided by: the recently updated Coupled Model Intercomparison Project Phase 5 (WORLDCLIM 1.4); testing data from Global Climate Models (GCM) CCSM4 (Community Climate System Model); MIROC-ESM (Model for Interdisciplinary Research on Climate); and the MPI-ESM-P (Max Planck Institute Earth system model) for the Mid-Holocene (30 second resolution, i.e. approximately 1 × 1 km) and LGM periods (2.5 min resolution, i.e. approximately 5 × 5 km; this layer was resampled to a 1 × 1 km resolution raster but not interpolated, meaning the resolution of interpretation for the LGM projected model and derived inferences must be 5 × 5 km). For the LIG chronozone, data from [START_REF] Otto-Bliesner | Simulating Arctic climate warmth and icefield retreat in the Last Interglaciation[END_REF] were used (30 second resolution).

Suitability values as modelled by MAXENT were converted into presence/absence, using the ten-percentile threshold of 0.305571, the minimum probability for presence cells after discarding the lowest 10% of the suitability distribution [START_REF] Pearson | Predicting species distributions from small numbers of occurrence records: a test case using cryptic geckos in Madagascar[END_REF]. Long-term persistence areas were inferred as continuously suitable areas (> ten-percentile value) over the LIG, LGM and Mid-Holocene periods.

PHYLOGEOGRAPHICAL RECONSTRUCTION OF MYRTUS COMMUNIS Genetic sampling strategy and molecular methods

We collected 118 samples from the entire Mediterranean and Macaronesian distribution area of M. communis. One randomly-sampled individual per site was newly genotyped using AFLP (Table S1; [START_REF] Vos | AFLP: a new technique for DNA fingerprinting[END_REF]. The AFLP reaction was performed on 300 to 500 ng of DNA extract, as described in [START_REF] Migliore | Genetic diversity and structure of a Mediterranean endemic plant in Corsica (Mercurialis corsica, Euphorbiaceae)[END_REF]. After screening selective primers, three primer combinations giving clearly visible band profiles were chosen (EcoRI-AAGG/MseI-CCAG, EcoRI-AAC/MseI-CAC, EcoRI-AAC/MseI-CAA). Polyacrylamide electrophoreses (0.4%) were performed using a 96-capillary automated sequencer (Megabace 1000, Amersham Bioscience), with manual scoring used to detect error-prone markers that had no reliable peak. The reproducibility of the AFLP markers was checked by repeating the complete analysis on 20 samples for each pair of primers [START_REF] Bonin | How to track and assess genotyping errors in population genetics studies[END_REF]. A genotyping error rate of 5.7% was obtained, due to the unreliability of some markers, which were subsequently removed from the dataset.

We also included sequence data for the plastid DNA trnL-trnF and rpl2-trnH intergenic regions, the nuclear DNA external transcribed spacer (ETS) region, and the internal transcribed spacer (ITS1-5.8S-ITS2) region, (n = 309 and 176, respectively; Table S1). The methods used to extract DNA and to sequence DNA markers are described in Migliore et al. (2012), since these sequences were previously generated for phylogenetic reconstruction and molecular dating. Genetic analyses AFLP genotypes were assigned to genetically homogeneous clusters via the model-based clustering algorithm provided by STRUCTURE 2.3.3 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF][START_REF] Falush | Inference of population structure using multilocus genotype data: dominant markers and null alleles[END_REF]. Bayesian analysis was run for 1,000,000 generations (burn-in of 100,000), and for K between 2 and 21 with ten iterations for each K-value. The admixture and recessive allele models were then chosen. The most likely number of clusters was determined using the criteria lnPr(X/K) and ΔK [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF].

After independently concatenating the two plastid and the two nuclear regions, plastid haplotype and nuclear ribotype median-joining networks were constructed, using NETWORK 4.613 [START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF].

To check the congruence of the genetic patterns from the same samples studied with AFLP and DNA sequences, a Mantel test was undertaken, comparing inter AFLP genotype distances (Jaccard index) to nucleotidic divergence of plastid DNA haplotypes and nuclear DNA ribotypes (p distance). The Mantel test was based on 10,000 permutations (GENALEX 6.5; [START_REF] Peakall | GENALEX 6.5: genetic analysis in Excel. Population genetic software for teaching and research-an update[END_REF].

Genetic diversity and singularity were summarized by several indices calculated using GENALEX 6.5: the number of plastid DNA haplotypes (NH) and nuclear DNA ribotypes (NR) including private haplotypes and ribotypes; the unbiased diversity (uh) for all the markers; the mean genetic distance within each group of samples sequenced (GD); the number of AFLP bands (NB); the number of private AFLP bands (PB); and the percentage of AFLP polymorphic loci (%P). All these indices were computed according to (i) the main AFLP genetic clusters detected, and (ii) the main geographical trends of climatic suitability found according to SDMs, excluding the Macaronesian genetic cluster whose sampling was too low.

RESULTS

MODELLED CURRENT AND PAST SUITABLE AREAS FOR MYRTUS COMMUNIS

The most parsimonious fitted model (based on AIC, Table S2, Figure 1B) was a combination of linear and quadratic features, with a True Skill Statistics (TSS) value of 0.73 (calculated using the 10-percentile threshold for a binary transformation; and an AUC value for training data of 0.88 ± 0.0068 (gain threshold reached after 1,040 iterations). Jackknife and response curves indicated that higher suitability values were inferred with high values of "minimal temperature of coldest month" with an optimum at 6°C (BIO6; test gain when used in isolation (Tg) of 0.88). Suitability also increased with moderate values of "mean temperature in the wettest quarter" with an optimum at 18°C (BIO8; Tg = 0.49) and low values of "annual temperature range" (optimum = 8°C for BIO7; Tg = 0.31). Annual precipitation values also positively influenced the suitability level (Tg = 0.25), with an inferred maximum suitability of 825 mm of annual precipitation. The variables "precipitation in the driest month" (BIO14) and "isothermality" (BIO3) slightly influenced the current model (Tg = 0.025 and 0.023 respectively). In summary, M. communis appeared to be particularly sensitive to cold stress, with suitable areas characterized by mean minimal temperature above frost threshold and mild winter temperatures. Areas with low level of precipitation (from 220 mm) are suitable for M. communis, but its optimum is around the maximum average annual precipitation value that defines the Mediterranean climate. The low number of available Myrtle fossils from the LGM and Mid-Holocene periods (Migliore et al., 2012) makes it difficult to choose between one of the three past climate models (CCSM4, MIROC-ESM, and MPI-ESM-P) for projecting SDM. The myrtle range dynamics over time using the MIROC-ESM climate model is the only one presented in this manuscript (Figure 2; results with the two other past climate models are presented as Supplementary Information). However, the SDM inference of long-term persistence areas is presented using all three past-climate model projections (Figure 3).

During the LIG period, M. communis was inferred to be potentially widespread throughout the entire Mediterranean Basin, albeit with a suitable area along the European and northern African coasts that was narrower than the current one (Figure 2A). The total inferred size of suitable areas approached 4,468,025 km² during the LIG versus 13,482,343 km² at present. A high level of suitability was suggested along the Atlantic coast (Portugal, Gibraltar, the south-western Iberian Peninsula, and western Morocco). In the eastern Mediterranean, the spatial range was similar to the current one, with a slight southward shift of suitable areas.

During the LGM, the whole study area was much less suitable for M. communis than it is under current climate conditions (Figure 2B). A drastic restriction of its potential suitable range is suggested: on average, the climate models showed 4,343,671 km 2 of suitable areas lost between the LIG and the LGM (Figures S1-2). We have to point out that the physiography of the Mediterranean Basin was quite different during these periods from how it is at present: in particular, the Adriatic Gulf had seen a sea-level drop of 100-120 m relative to the present. Using the LIG map as a baseline for analysing the changes occurring during the LGM, we can infer that the suitable areas for M. communis could have moved down towards lower altitudes and southward towards lower latitudes. This was apparent on western and central Mediterranean islands, along the Atlantic coastline, and, to a lesser degree, in Greek and Turkish areas (Figure 2B). Among the climate models used for the LGM, the CCSM4 model was the most restrictive with only 91,052 km² of suitable area, in contrast to 120,428 km² for the MIROC-ESM (with higher suitability in Iberian Peninsula and the Gibraltar Strait), and 161,585 km² for the MPI-ESM-P models (with higher suitability in western Mediterranean islands) (Figures S1-2).

During the Mid-Holocene, the climatic conditions were inferred to be more suitable than during the LGM, thus allowing for a M. communis range expansion (Figures 2C,S3,S4). Averaging over the three climate models, this gain in area comes to 196,412 km² in total.

An important but unexpected result was the low number and small size of inferred longterm refugia (i.e. continuously suitable areas over the LIG, LGM and Mid-Holocene), whatever past-climate model is used (Figure 3). We infer drastic changes in the spatial range of the species between these periods together with a quasi-absence of overlap between suitable areas during the three contrasted climatic periods, apart from parts of the Macaronesian islands (Azores and Madeira), around the Gibraltar Strait, on the island of Crete, and in the Levant (Figure 3). Above the 41°N parallel, no long-term persistence areas over the LIG, LGM and Mid-Holocene were detected in the northern part of the Mediterranean Basin (Figure 3).

MULTI-MARKERS GENETIC DIVERSITY OF MYRTUS COMMUNIS

The AFLP analysis on 118 individuals of M. communis resulted in 199 scored fragments, 191 of which were found to be polymorphic, with a length varying from 66 to 460 bp (dataset available on request). The most probable partitioning of the genetic variation of M. communis according to STRUCTURE was in K = 2 and K = 4 genetic clusters (Figures 4A-B, S5).

The genetic distances calculated from AFLP, plastid and nuclear DNA sets of data were positively correlated and the Mantel tests were significant for the Mediterranean samples (n = 110: plastid DNA vs AFLP: R² = 0.20, P = 0.0001; plastid DNA vs nuclear DNA: R² = 0.19, P = 0.0001; nuclear DNA vs AFLP: R² = 0.20, P = 0.0001). Genetic diversity indices were computed for each AFLP genetic cluster: for the three molecular markers, the Western Mediterranean (WM) and Eastern Mediterranean (EM) clusters had similar levels of diversity when accounting for the different levels of sampling (uh; Table 1). The circum-Mediterranean (CM) cluster had the lowest diversity of plastid DNA markers but had the highest diversity of nuclear DNA, and it also had the highest percentage of polymorphic markers for AFLP (%P; Table 1). This contrast could be due to the intermediate geographical position of the CM cluster providing a higher admixture of the WM or EM clusters, i.e. contact zones (see below).

The plastid DNA network was characterized by 12 haplotypes (Figures 4A-C, S6A-C) geographically structured within the Mediterranean Basin and between the Macaronesian and Mediterranean regions. The simultaneous analysis of ETS and ITS nuclear DNA data revealed a similar genetic pattern to that obtained with plastid DNA data; 41 ribotypes were distributed through the three main genetic clusters (Figure 4A-D and Figure S6B-D).

Higher values of genetic singularity and diversity were detected south of the 41°N parallel (Figure 3), as follows: 9 plastid DNA haplotypes (vs 4 to the north) of which 5 (vs 0) are private; 26 ribotypes (vs 19) of which 17 (vs 11) are private; 175 AFLP bands (vs 135) of which 50 (vs 10) are private (Table 1). Higher genetic diversity in the southern Mediterranean was detected similarly for plastid DNA, nuclear DNA and AFLP molecular markers (see uh and %P in Table 1). Southern France and Corsica, which were the most-sampled areas in the northern Mediterranean, appeared to have a low genetic originality (no private markers) even though France was a contact zone between the WM and CM genetic clusters (Figures 4 andS6). Considering each of these clusters independently, we found that the areas harbouring the highest number of plastid DNA haplotype were situated in North Africa (WM, Algeria and Tunisia), in Sicily (CM), and in the Levant and Cyprus (EM). The number of ribotypes was also higher in North Africa (WM), the Levant, and Cyprus (EM) (Figure S6).

At the scale of the Mediterranean Basin, the contact zones between genetic clusters occurred longitudinally, and originated from WM and CM clusters, whereas plastid haplotypes and ribotypes of the EM cluster remained isolated (Figure 4). The major plastid DNA contact zone was detected between EM and CM clusters in Turkey and the Levant, implying widespread haplotypes 4 and 12 (Figure S6). The same pattern was observed for nuclear DNA markers (AFLP and sequences), except that the WM cluster was detected in Turkey and the Levant. In summary, the three molecular markers are consistent in indicating that the Western and Central parts of the Mediterranean could have been stronger sources of range expansions than the Eastern parts.

DISCUSSION

Palaeoclimatic SDM and levels of genetic diversity suggest drastic changes in suitable areas for M. communis during the last 130,000 years. In contrast to the phylogeography of Olea europaea (Besnard et al., 2013), past suitability modelling hindcasts the persistence of myrtle populations in a few small long-term refugia.

Projected past suitable areas for M. communis indicated a strong signature of the last glaciation on the potential geographical range of the species. On average 97% of suitable areas were inferred to be lost between the LIG and the LGM (Figure 2). The decrease in genetic diversity and singularity revealed by molecular markers (Table 1) is consistent with the idea that northern populations along the Mediterranean coastline were the most affected by LGM climate change, as also evidenced by Laurus nobilis (Rodríguez-Sánchez & Arroyo, 2008), Erica scoparia (Désamoré et al., 2012), and Chamaerops humilis (García-Castaño et al., 2014). Furthermore, dramatic consequences of glaciations with major bottlenecks, or even complete extinction in the Mediterranean, were suggested for some thermophilous plants, as Nerium oleander L. (Mateu-Andrés et al., 2015), and Ceratonia siliqua L. [START_REF] Ramón-Laca | The ecological status of the carob-tree (Ceratonia siliqua, Leguminosae) in the Mediterranean[END_REF].

The strong range contraction inferred for M. communis during the LGM raises the issue of its long-term survival through the Pleistocene. Since fossil data are often too scarce for formally validating SDM, the identification of refugia using SDM approaches remains tentative [START_REF] Ashcroft | Identifying refugia from climate change[END_REF]. According to SDM, the size of the myrtle's long-term refugia is highly variable with (i) two main poles of persistence in the western Mediterranean (Atlantic coast and Baetic-Rifan complex) and in the east (Levantine coast), and (ii) small pockets along the mainland and island coastlines (Figure 3). The reduced overlap of suitable areas among the three palaeomaps (LIG, LGM, and Mid-Holocene) leads to a few very reduced long-term refugia compared to the size of the glacial refugia (Figures 23). These results, rarely found in the literature, question the long-term refugia hypothesis and suggest the possibility of a regional persistence of M. communis by range shifts toward temporary refugia, where conditions became suitable. It must be noted that the projections of SDM at the LIG, LGM and Mid-Holocene periods do not aim to infer M. communis distribution at those times, but rather to identify potential suitable areas considering the environmental variables selected. Distribution and suitable areas would be consistent only under the assumption of unlimited dispersal of the species.

At this juncture, we may examine our capacity to reject the possibility of long-term and in situ refugia for M. communis that could be cryptic for our methods. The role of local microrefugia may have played a more important role than expected, as the heterogeneity of micro-habitats may have allowed species to cope with environmental stress at finer spatial scales [START_REF] Stewart | Cryptic northern refugia and the origins of the modern biota[END_REF][START_REF] Serra-Diaz | Disturbance and climate microrefugia mediate tree range shifts during climate change[END_REF][START_REF] Meineri | Fine-grain, large-domain climate models based on climate station and comprehensive topographic information improve microrefugia detection[END_REF]. However, it is worth considering the studies conducted on M. communis at small scales. Studies in southern Spain are consistent with a low likelihood of long-term persistence of M. communis within small remnant patches of woodlands [START_REF] González-Varo | Extinction debt of a common shrub in a fragmented landscape[END_REF]. The regeneration and genetic diversity of M. communis are generally lower in small populations, and have a high sensitivity to fragmentation (e.g. [START_REF] González-Varo | Effects of fragmentation on pollinator assemblage, pollen limitation and seed production of Mediterranean myrtle (Myrtus communis)[END_REF][START_REF] Nora | Genetic variation and structure in the Mediterranean shrubs Myrtus communis and Pistacia lentiscus in different landscape contexts[END_REF]. Though places such as Corsica can harbour conditions for local microrefugia [START_REF] Médail | Glacial refugia influence plant diversity patterns in the Mediterranean Basin[END_REF][START_REF] Gavin | Climate refugia: joint inference from fossil records, species distribution models and phylogeography[END_REF], local persistence is unlikely, perhaps due to the severe effects of glaciations [START_REF] Conchon | Quaternary glaciations in Corsica[END_REF]. Despite its abundance in the lowlands of Corsica up to 500 m of altitude, and irrespective of which molecular marker was used, we found a lack of genetic differentiation, with for example only two plastid DNA haplotypes (8, 9) detected after an extensive sampling. These haplotypes are not private to Corsica and belong to the WM genetic cluster, occurring in Sardinia, South-east France, Italy and Eastern Algeria (Figures 4 andS6). Moreover, all the 12 Corsican palaeoecological records for Myrtus correspond to the late Holocene [START_REF] Reille | Quelques aspects de l'activité humaine en Corse durant le Subatlantique et ses conséquences sur la végétation[END_REF][START_REF] Reille | New pollen-analytical researches in Corsica: the problem of Quercus ilex L. and Erica arborea L., the origin of Pinus halepensis Miller forests[END_REF][START_REF] Reille | New pollen-analytical researches in Corsica: the problem of Quercus ilex L. and Erica arborea L., the origin of Pinus halepensis Miller forests[END_REF][START_REF] Vella | Approches géomorphologique et géophysique des interactions sociétés/milieux en Corse au cours de l'Holocène[END_REF]. [START_REF] Reille | New pollen-analytical researches in Corsica: the problem of Quercus ilex L. and Erica arborea L., the origin of Pinus halepensis Miller forests[END_REF] suggested that M. communis is "the latest arrival in the thermophilous vegetation directly related to human action".

There are caveats to our questioning of long-term in situ refugia for M. communis. Specifically, reports that fleshy-fruited myrtle berries are mainly dispersed by frugivorous birds [START_REF] Herrera | A study of avian frugivores, bird-dispersed plants, and their interaction in Mediterranean scrublands[END_REF][START_REF] Traveset | Ecology of fruit-colour polymorphism in Myrtus communis and differential effects of birds and mammals on seed germination and seedling growth[END_REF][START_REF] Gonzalez-Varo | Fragmentation, habitat composition and the dispersal/predation balance in interactions between the Mediterranean myrtle and avian frugivores[END_REF] are consistent with an active role for migration. Moreover, widespread plastid DNA haplotypes and ribotypes (Figure S6) suggest episodes of frequent and/or long-distance migrations of M. communis in the western and central parts of the Mediterranean Basin. Our analyses revealed that migrations have been less pronounced in the eastern Mediterranean. Genetic data also suggest that migrations occurred rapidly, recently, and/or over long distances to prevent molecular divergence between remote sites (e.g. the CM genetic cluster).

The phylogeography of Smilax aspera L., a climbing plant with sclerophyllous leaves and red berries dispersed by birds, is also congruent with the pattern described for M. communis with a longitudinal rather than a latitudinal structure of genetic diversity [START_REF] Chen | Understanding the formation of Mediterranean-African-Asian disjunctions: evidence for Miocene climate-driven vicariance and recent long-distance dispersal in the Tertiary relict Smilax aspera (Smilacaceae)[END_REF]. Such a longitudinal organization of phylogroups more or less isolated, but rarely totally disjunct constitutes thus one major phylogeographical trend of the circum-Mediterranean woody species (see the review of [START_REF] Feliner | Patterns and processes in plant phylogeography in the Mediterranean Basin. A review[END_REF]. As for M. communis, the phylogeography of Erica arborea (Désamoré et al., 2011), Arbutus unedo (Santiso et al., 2016), and Laurus nobilis [START_REF] Rodriguez-Sanchez | Late Neogene history of the laurel tree (Laurus L., Lauraceae) based on phylogeographical analyses of Mediterranean and Macaronesian populations[END_REF]) support a stronger diversification within the western Mediterranean, and inferred eastward migrations for Myrtus and Arbutus (see Santiso et al., 2016).

CONCLUSION

Documentation of past range dynamics of Mediterranean woody taxa supports the idea that glaciations led to significant changes in Mediterranean vegetation, such as in the abundance of most thermophilous plants. Our results are in accordance with recent findings indicating that the Mediterranean Basin was less stable from a palaeoenvironmental point of view, having experienced higher extinction rates of species and lineages compared to other Mediterraneantype-ecosystems [START_REF] Valente | Contrasting evolutionary hypotheses between two mediterraneanclimate floristic hotspots: the Cape of southern Africa and the Mediterranean Basin[END_REF][START_REF] Cowling | Variation in plant diversity in mediterraneanclimate ecosystems: the role of climatic and topographical stability[END_REF]. Myrtus communis is representative of a Tertiary cold-sensitive lineage and we propose here that the severe effects of Pleistocene glaciations were faced despite a relatively reduced role of long-term refugia. Myrtus communis seems to have survived the Quaternary glaciations by regional range shifts towards temporary refugia, a response which has, to our knowledge, not been considered in previous Mediterranean phylogeographical or SDM studies. This hypothesis is an important issue regarding conservation or restoration of connectivity between Myrtus populations, and more generally between fragments of woodlands to ensure the future resilience of Mediterranean thermophilous shrublands.

Table S3. Bibliographic review of the main imprints left by past environmental changes on Mediterranean widespread, Tertiary-originated, and/or thermophilous woody species, combining genetic, distribution modelling and palaeoecological data.

Bibliographic review of the main imprints left by past environmental changes on Mediterranean widespread, Tertiary-originated, and/or thermophilous woody species, combining genetic, distribution modelling and palaeoecological data. Mean values of cold resistance for vegetative organs in adult specimens are those causing 50% damage to leaf buds and cambium (based on [START_REF] Larcher | Low temperature effects on Mediterranean sclerophylls: an unconventional viewpoint[END_REF][START_REF] Larcher | Temperature stress and survival ability of Mediterranean sclerophyllous plants[END_REF][START_REF] Quézel | Ecologie et biogéographie des forêts du bassin méditerranéen[END_REF][START_REF] Flexas | Photosynthetic limitations in Mediterranean plants: a review[END_REF]. Potential distribution of Myrtus communis fitted to current climatic conditions from actual occurrences using MAXENT and projected during the Last Inter-Glacial period (A), and the Last Glacial Maximum using the CCSM4 (B), MIROC-ESM (C) and MPI-ESM-P (D) climate models. Continuous suitability values modelled with MAXENT were converted into presence/absence, using the ten percentile threshold method (see Methods). Shades referred to palaeodistributions modelled just before the chronozone represented. Table S2. AIC scores and number of parameters for fitted MAXENT models.

Akaike Information Criterion (AIC) scores and number of parameters for all fitted MAXENT models using all possible simple and combination of features except "threshold": L = linear, Q = quadratic, P = product. Table S3. Bibliographic review of the main imprints left by past environmental changes on Mediterranean widespread, Tertiary-originated, and/or thermophilous woody species, combining genetic, distribution modelling and palaeoecological data.

Features combination

Bibliographic review of the main imprints left by past environmental changes on Mediterranean widespread, Tertiary-originated, and/or thermophilous woody species, combining genetic, distribution modelling and palaeoecological data. Mean values of cold resistance for vegetative organs in adult specimens are those causing 50% damage to leaf buds and cambium (based on [START_REF] Larcher | Low temperature effects on Mediterranean sclerophylls: an unconventional viewpoint[END_REF][START_REF] Larcher | Temperature stress and survival ability of Mediterranean sclerophyllous plants[END_REF][START_REF] Quézel | Ecologie et biogéographie des forêts du bassin méditerranéen[END_REF][START_REF] Flexas | Photosynthetic limitations in Mediterranean plants: a review[END_REF]. -Evidence for Pleistocene expansion and admixture.

-Range contraction stronger for LIG than LGM, mainly for the continental areas of France-Italy and E Spain-Sardinia-NE Morocco.

-Very low likelihood of occurrence for Corsica contrary to other islands.

-Leaf buds: -11.5°C -Two waves of expansions in the Miocene and Pleistocene.

-One widespread haplotype supports recent events of long distance dispersal.

-Migrations to Macaronesia and Sahara (Tibesti Mt).

-Range contraction during the LGM.

-No suitable area predicted for LGM in W Europe but refugia along the African Atlantic coast.

- Désamoré et al., 2011Désamoré et al., , 2012 Laurus nobilis L.

( -The WM genetic cluster is more diverse and expanded eastward.

-One widespread haplotype supports recent events of long distance dispersal.

-Large suitable areas during the Pliocene -EM genetic cluster: higher genetic diversity and stronger genetic structure.

-Bottleneck during the LGM followed by rapid post-glacial recolonisations suggested to explain the lower diversity and differentiation in the WM. al., 2007[START_REF] Magri | The distribution of Quercus suber chloroplast haplotypes matches the palaeogeographical history of the western Mediterranean[END_REF][START_REF] Vessella | Quercus suber range dynamics by ecological niche modelling: from the Last Interglacial to present time[END_REF]

Figure 1 .

 1 Figure 1. Known occurrences of Myrtus communis in the Mediterranean Basin, gathered from botanical databases and literature in addition to field collected data (n = 4,168) (A). Current suitable areas for the species, applying species model distribution (SDM) from actual occurrences using MAXENT with six bioclimatic variables (see Methods) (B).

Figure 2 .

 2 Figure 2. Suitable areas for Myrtus communis (projection of the most parsimonious MaxEnt model fitted to current climatic conditions from actual occurrences) during the Last Inter-Glacial (LIG) period (A), the Last Glacial Maximum (LGM) using the MIROC-ESM climate model (B), the Mid-Holocene using the MIROC-ESM climate model (C), and the Present (D). Continuous suitability values modelled with MAXENT were converted into presence/absence, using the ten-percentile threshold ≥ 0.3 (see Methods and Figures S1-4). Shading refers to palaeodistributions modelled just before the chronozone represented.

Figure 3 .

 3 Figure 3. Long-term persistence areas inferred as continuously suitable areas over the Last Inter-Glacial (LIG) period, the Last Glacial Maximum (LGM) and Mid-Holocene using the CCSM4 (A), MIROC-ESM (B) and MPI-ESM-P (C) climate models (ten-percentile threshold ≥ 0.3; 5 × 5 km resolution).

Figure 4 .

 4 Figure 4. Spatial genetic structure of Myrtus communis through its Mediterranean distribution range (A), based on the AFLP (Amplified Fragment Length Polymorphism) data summarized from STRUCTURE, based on K = 4 genetic clusters (B), and the median-joining networks from concatenated plastid DNA sequences trnL-trnF and rpl2-trnH intergenic regions (C) and concatenated nuclear DNA sequences external transcribed spacer (ETS) region and internal transcribed spacer (ITS1-5.8S-ITS2) (D). Colour chart is based on the four main genetic clusters identified by AFLP: EM (Eastern Mediterranean in orange), WM (Western Mediterranean in green), and CM (circum-Mediterranean in blue), and Macaronesia (in purple). Numbers in the panel (A) referred to spatially-closed populations whose samples have the same genetic clusters identified. Haplotype and ribotype identifiers are detailed in Figure S6; circle sizes in networks are proportional to haplotype/ribotype frequencies.

Figure S1 .

 S1 Figure S1. Palaeodistribution of Myrtus communis during the LIG and LGM. MAXENT projection of the model of current distribution of Myrtus communis during the Last Inter-Glacial period (ca. 116-130 ka BP) (A), and the Last Glacial Maximum (ca. 19-26.5 ka BP) according to the CCSM4 (B), MIROC-ESM (C) and MPI-ESM-P (D) climate models.

Figure S2 .

 S2 Figure S2. Presence/absence of Myrtus communis during the LIG and LGM.

Figure S3 .

 S3 Figure S3. Palaeodistribution of Myrtus communis during the Mid-Holocene. MAXENT projection of the model of current distribution of Myrtus communis during the Mid-Holocene (ca. 6 ka BP) according to the CCSM4 (A), MIROC-ESM (B) and MPI-ESM-P (C) climate models, and during the Present (D).

Figure S4 .

 S4 Figure S4. Presence/absence of Myrtus communis during the Mid-Holocene.Potential distribution of Myrtus communis fitted to current climatic conditions from actual occurrences using MAXENT and projected during the Mid-Holocene using the CCSM4 (A), MIROC-ESM (B) and MPI-ESM-P (C) climate models, and during the Present (D). Continuous suitability values modelled with MAXENT were converted into presence/absence, using the ten percentile threshold method (see Methods). Shades referred to palaeodistributions modelled just before the chronozone represented.

Figure S5 .

 S5 Figure S5. Analysis of STRUCTURE results for AFLP data.Representation of the Log probability of data L(K) (A) and the magnitude of ΔK (B) as a function of K for the AFLP dataset of Myrtus communis, according to[START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF] and the software STRUCTUREHARVESTER(Earl et al., 2012).

Figure S6 .

 S6 Figure S6. Spatial genetic structure with haplotypes identifiers of Myrtus communis, based on plastid DNA and nuclear DNA. Spatial genetic structure with haplotypes identifiers of Myrtus communis through its Mediterranean distribution range, based on the median-joining networks from plastid DNA sequences trnL-trnF and rpl2-trnH intergenic regions (A,C) and nuclear DNA sequences external transcribed spacer (ETS) region and internal transcribed spacer (ITS1-5.8S-ITS2) (B,D). When more than one individual was sequenced per population circles on the maps are horizontally divided in two, and if spatially-closed populations have the same haplotypes or ribotypes they are separated by the symbol "/". Circle sizes in networks are proportional to haplotype/ribotype frequencies.

  (2,401-2,449 bp: trnH(GUG)-psbA, trnL(UAG)-rpl32, trnT(UGU)-trnL(UAA) and Intron L; n = 207, 23 pops) > 28 haplotypes -Molecular dating -(Map of LGM climatic conditions) 2 main genetic clusters -Tertiary divergence between A. unedo (WM) / A. andrachne (EM) -Western diversification of A. unedo during the last 700 ka (hardest Pleistocene glaciations) and eastward migrations. -Key role of glacial refugia (WM + NAfrica) + Refugia within refugia model. -Recent post-glacial migration to Ireland from Iberia662 bp: ITS; n = 15, 15 pops with cloning) > 19 ribotypes -Poor phylogeographical resolution of nrDNA ITS sequences. EM (Anatolia) retained ancestral ribosomal sequences. in the WM, potentially related to a Miocene W/E vicariance (5.83-8.32 Ma).

-

  Pliocene divergence of 3 genetic clusters (5.82-1.49 Ma) with a E/W bipolarisation of genetic diversity.

-

  ) -2 pDNA (2,562 bp: trnK-MatK and trnD-T; n = 57) Longitudinal geographical structure related to isolation of WM and EM genetic clusters with strong barriers in the East (Floristic Rechinger's line).

-

  (932 bp: trnL-trnF and rpl2-trnH; n = 309) > 12 haplotypes -2 nrDNA (1,011 bp: ITS and ETS; n = 176) >39 ribotypes -2 AFLP (n = 118) -Molecular dating -CMIP5 Worldclim, MaxEnt (GPS n = 4,168): LIG, LGM, Mid-Holocene, Present -Pleistocene pollen (n = 76) -Neogene fossils 3 main clusters -A vicariance event in the late Miocene led to W/E bipolarisation of genetic diversity. -Higher diversification in WM with eastward migration. CM genetic cluster with potentially a central origin, migrated westward and eastward (one widespread haplotype). EM genetic cluster remained geographically restricted. -Migrations to Macaronesia and Sahara (Myrtus nivellei). supposed to explain the absence of pDNA variability. -Unique divergence detected Sahara/ Med: 7.2-1.2 Ma: migrations from the Saharan refugia? -High sensitivity of riparian corridors to frost during the LGM supposed. Divergence between Oleaster lineages started during the late Pliocene. -Despite long distance dispersal and admixture with cultivated Olives, a strong W/E differentiation is still present. Several haplotypes restricted to the East. One genetic cluster could have central origin. -Migrations/introgression in Macaronesia/Sahara. RFLP (n = 570, 57 pops) > 3 mitotypes -5 cpSSRs (n = 1339, 48 pops) and genetic interpolation -Worldclim, GLMs /GAMs /RF /CTA /MaxEnt (GPS n = 128,653): Present 3 main genetic clusters -Eight gene pools organized in three main genetic clusters. -Post glacial expansions from a minimum of 3 refugia to explain the observed phylogeography and pattern of population genetic diversity. SSR (n = 856 , 34 pops) No pDNA variability -Demographic bottleneck suggested: 1 WM haplotype and 2 in the EM.-After the bottleneck long distance dispersal or human-mediated dispersal could explain the current Mediterranean range. organized longitudinally, with an eastern one, centred on Anotalia, inferred as ancestral.-The split between the Eastern (Anotalia) and Central (Balkan) genetic clusters occurred before in the Early Pliocene whereas the split between the central and western (Italy) genetic clusters is dated to the Mid Pleistocene. Not available -Bagnoli et al., 2016 50 -Haplotype sharing between Q. suber and Q. cerris in Pleistocene Italian glacial refugia, as hypothesis. (594 bp: ITS; n = 6) -5 AFLP (n = 19, 19 pops) -7 allozymes (n = 609, 24 pops) 2 main genetic clusters -Genetic clusters are corresponding to Q. calliprinos in EM / Q. coccifera in WM.

  (592 bp: ITS; n = 11) -5 AFLP (n = 41, 39 pops) 2 main genetic clusters -Strong imprint of glaciations with two main genetic clusters corresponding to the main peninsula refugia. -Admixture and post-glacial contact zone of lineages in the Rhône Valley (France). Extensive hybridization with Q. coccifera and Q. suber. 5 pDNA PCR RFLP (n ≈ 365, 91 pops/n = 854, 94 pops) -1 nrDNA (617 bp: ITS; n = 17) -5 AFLP (n = 34, 32 pops) -PMIP3 Worldclim, openModeller algorithms (GPS n = 63,733): LIG, LGM, Mid-Holocene, Present -Pleistocene pollen and tree fossils (n = 17) 3 main genetic clusters -3 genetic clusters (W. Iberia, Morocco / Algeria, Tunisia, Tyrrhenian islands / Italy, Sicily) considered under a hypothesis of deep divergence that could be related to Oligocene-Miocene palaeogeography, but a molecular dating approaches is lacking to test it. -Long-term persistence in Mediterranean without striking range contraction or shift, even during the LGM. -Leaf buds: -16°C -Cambium: -26°C -Lumaret et al., 2005 -Lopez de Heredia et

Table 1 .

 1 Genetic (plastid DNA sequences, nuclear DNA sequences and AFLP data) features for Mediterranean samples of Myrtus communis in relation with the three detected genetic clusters (Eastern, Circum and Western Mediterranean) and with the distribution of samples in the North or in the South of the 41°N parallel: sample size (n), number of plastid DNA haplotypes (NH with number of private haplotypes in brackets), unbiased diversity (uh with standard errors in brackets), mean genetic distance (GD), number of nuclear DNA ribotypes (NR with number of private nrDNA ribotypes in brackets), number of AFLP bands (NB), number of private AFLP bands (PB), percentage of polymorphic loci (%P).

		Plastid DNA sequences	Nuclear DNA sequences			AFLP data
		n	NH	uh	GD	n	NR	uh	GD	n	NB	PB	uh	%P
	Cluster EM Eastern Mediterranean	60	2	0.051 (0.051)	0.463	26	6	0.042 (0.021)	6.108	10	110	12	0.145 (0.013)	43.72
	Cluster CM Circum Mediterranean	86	3	0.037 (0.028)	0.333	52	16	0.091 (0.030)	7.177	48	151	23	0.124 (0.011)	67.34
	Cluster WM Western Mediterranean	155	4	0.055 (0.032)	0.491	94	15	0.041 (0.018)	3.831	55	132	19	0.130 (0.013)	54.27
	All data	301	9	0.212 (0.057)	/	172	37	0.294 (0.035)	/	113	199	/	0.278 (0.018)	/
	NM Northern Mediterranean (41°N)	133	4 (0)	0.085 (0.048)	0.768	86	19 (11)	0.186 (0.032)	8.332	43	135	10	0.161 (0.014)	56.78
	SM Southern Mediterranean (41°N)	168	9 (5)	0.267 (0.061)	2.407	86	26 (17)	0.340 (0.036)	11.844	70	175	50	0.173 (0.013)	82.41

Table S1 . List of samples of Myrtus communis used for genetic analyses (plastid and nuclear sequencing, and AFLP genotyping).

 S1

	Samples Locations ALBA1 Albania Morocco MRIB5 Morocco PALE1 Portugal PALE9 Portugal ASMA1 Portugal (Azores) ASMA9 Portugal (Azores) ALBA9 MRIB1 ASMP1 Portugal (Azores)	X 19.483458 -5.69763000 -5.69763000 -7.65744 -7.65744 -25.161912 -25.161912 -25.105245	Y 40.98658 34.79624500 34.79624500 37.653408 37.653408 36.967147 36.967147 36.953679	AFLP clusters CM / WM WM / Mac / Mac	Plastid DNA haplotypes 4 14 14 14 14 6 7 7	GenBank trnL-trnF FJ611337 HM014154 GenBank rpl2-trnH FJ611338 HM014150 FJ611338 HM014150 FJ611338 HM014150 FJ611338 HM014150 FJ611338 HM014147 FJ611338 HM014146 FJ611338 HM014146	Nuclear DNA ribotypes 16 26 / 31 / 38 / 39	GenBank ITS1-5,8S-ITS2 GU984346 HM106407 GenBank ETS-18S GU984342 HM106410 / / GU984340 HM106410 / / JF304834 HM106419 / / GU984370 HM106421
	ASMP9	Portugal (Azores)	-25.105245	36.953679	/	9	FJ611338 HM014149	/	/	/
	MABO1	Portugal (Madeira)	-16.771892	32.752919	Mac	13	FJ611338 HM014153	40	GU984370 HM106417
	MABO5	Portugal (Madeira)	-16.771892	32.752919	/	13	FJ611338 HM014153	/	/	/
	MAFA1	Portugal (Madeira)	-16.851667	32.78445	Mac	/	/	/	/	/	/
	MAVI0	Portugal (Madeira)	-17.023961	32.731635	Mac	13	FJ611338 HM014153	41	JF304831	JF304872
	ECAD1	Spain	3.28451100	42.30609200	/	9	FJ611338 HM014149	31	GU984340 HM106410
	ECAD5	Spain	3.28451100	42.30609200	WM	9	FJ611338 HM014149	/	/	/
	ECRE1	Spain	3.29008333	42.31738889	/	9	FJ611338 HM014149	26	JF304841 HM106410
	ECRE5	Spain	3.29008333	42.31738889	/	9	FJ611338 HM014149	/	/	/
	EDEV1	Spain	-0.31230556	39.34511111	/	8	FJ611338 HM014153	16	JF304853	JF304884
	EDEV5	Spain	-0.31230556	39.34511111	CM	4	FJ611337 HM014154	/	/	/
	EGUA10	Spain	-6.68666700	37.07972200	/	14	FJ611338 HM014150	/	/	/
	EGUA5	Spain	-6.68666700	37.07972200	WM	14	FJ611338 HM014150	31	GU984340 HM106410
	EMUR1	Spain	-0.35805556	39.12916667	/	8	FJ611338 HM014153	26	GU984342 HM106410
	EMUR5	Spain	-0.35805556	39.12916667	CM	8	FJ611338 HM014153	/	/	/
	EROS1	Spain	3.20763889	42.24469444	/	9	FJ611338 HM014149	30	JF304838 HM106410
	EROS5	Spain	3.20763889	42.24469444	CM	4	FJ611337 HM014154	/	/	/
	ESAL1	Spain	-0.32838889	39.38594444	/	4	FJ611337 HM014154	7	JF304819	JF304890
	ESAL5	Spain	-0.32838889	39.38594444	CM	14	FJ611338 HM014150	/	/	/
	ETOU1	Spain	-0.66052778	39.16016667	/	8	FJ611338 HM014153	20	JF304844	JF304884
	ETOU5	Spain	-0.66052778	39.16016667	CM	8	FJ611338 HM014153	/	/	/
	EESP1	Spain (Majorca)	3.37327831	39.70425638	/	8	FJ611338 HM014153	10	GU984353 JF304886

-Rodríguez-Sánchez & Arroyo, 2008 -Rodríguez-Sánchez et al., 2009

-[START_REF] Besnard | On the origin of the invasive olives (Olea europaea L., Oleaceae)[END_REF][START_REF] Nogués-Bravo | Predicting the past distribution of species climatic niches[END_REF], 2013[START_REF] Díez | Olive domestication and diversification in the Mediterranean Basin[END_REF][START_REF] Besnard | Single vs multiple independent olive domestications: the jury is (still) out[END_REF][START_REF] Díez | The jury may be out, but it is important that it deliberates: a response to Besnard and Rubio de Casas about olive domestication[END_REF]

-[START_REF] Burban | Phylogeography of maritime pine inferred with organelle markers having contrasted inheritance[END_REF] -Bucci et al., 2007 -Serra-Varela et al.,

ACKNOWLEDGEMENTS

This work was supported by the Provence-Alpes-Côte d'Azur Region, the Conservatoire Botanique National de Corse, and the Conservatoire Botanique National Méditerranéen de Porquerolles. Molecular biology analyses were carried out at the Molecular Biology Common Service of IMBE (Arbois) and the Molecular Biology Platform of INRA URFM (Avignon), with special thanks to Bruno Fady and Anne Roig. We are grateful to Michelle Barbier-Leydet for her help with the European Pollen Database. The English text was edited by Kolo Wamba (Proediting, Redwood City, USA). We thank the numerous contributors to the Myrtus sampling, cited in Migliore et al. (2012). Finally, we thank the Associate Editor Božo Frajman, Gonzalo Nieto Feliner and another reviewer for their very constructive comments, which helped to improve this manuscript.

SUPPORTING INFORMATION BOTANICAL JOURNAL OF THE LINNEAN SOCIETY

Surviving to glaciations in the Mediterranean region:

an alternative to the long-term refugia hypothesis *Corresponding author. E-mail: jeremy.migliore@imbe.fr / migliore.jeremy@libertysurf.fr

The following Supporting Information is available for this article: MAXENT projection of the model of current distribution of Myrtus communis during the Last Inter-Glacial period (ca. 116-130 ka BP) (A), and the Last Glacial ka BP) according to the CCSM4 (B), MIROC-ESM (C) and MPI-ESM-P (D) climate models.

Figure S2. Presence/absence of Myrtus communis during the LIG and LGM.

Potential distribution of Myrtus communis fitted to current climatic conditions from actual occurrences using MAXENT and projected during the Last Inter-Glacial period (A), and the Last Glacial Maximum using the CCSM4 (B), MIROC-ESM (C) and MPI-ESM-P (D) climate models. Continuous suitability values modelled with MAXENT were converted into presence/absence, using the ten percentile threshold method (see Methods). Shades referred to palaeodistributions modelled just before the chronozone represented. Potential distribution of Myrtus communis fitted to current climatic conditions from actual occurrences using MAXENT and projected during the Mid-Holocene using the CCSM4 (A), MIROC-ESM (B) and MPI-ESM-P (C) climate models, and during the Present (D). Continuous suitability values modelled with MAXENT were converted into presence/absence, using the ten percentile threshold method (see Methods). Shades referred to palaeodistributions modelled just before the chronozone represented.

Figure S5. Analysis of STRUCTURE results for AFLP data.

Representation of the Log probability of data L(K) (A) and the magnitude of ΔK (B) as a function of K for the AFLP dataset of Myrtus communis, according to [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF] and the software STRUCTUREHARVESTER (Earl et al., 2012). Spatial genetic structure with haplotypes identifiers of Myrtus communis through its Mediterranean distribution range, based on the median-joining networks from plastid DNA sequences trnL-trnF and rpl2-trnH intergenic regions (A,C) and nuclear DNA sequences external transcribed spacer (ETS) region and internal transcribed spacer (ITS1-5.8S-ITS2) (B,D). When more than one individual was sequenced per population circles on the maps are horizontally divided in two, and if spatially-closed populations have the same haplotypes or ribotypes they are separated by the symbol "/". Circle sizes in networks are proportional to haplotype/ribotype frequencies.

Table S1. List of samples of Myrtus communis used for genetic analyses (plastid and nuclear sequencing, and AFLP genotyping).

Table S2. AIC scores and number of parameters for fitted MAXENT models.

Akaike Information Criterion (AIC) scores and number of parameters for all fitted MAXENT models using all possible simple and combination of features except "threshold": L = linear, Q = quadratic, P = product.


Yann Goueffon
email: goueffon@chimie.ups-tlse.fr

Muguette Monjauze

Catherine Mabru

Laurent Arurault

Claire Tonon

Pascale Guigue

STUDY OF DEGRADATION MECHANISMS OF BLACK ANODIC FILMS IN SIMULATED SPACE ENVIRONMENT

Black anodic films are used on aluminium alloys in space applications to provide specific thermo-optical properties to the surface. During the process, thermal and/or residual stresses can lead to crazing of the film in particular conditions. After thermal cycling performed to simulate the space environment, cases of film flaking have been observed. This constitutes potential risks for the satellite lifetime. Peel-test and scratch-test have been used to evaluate the evolution of the adhesion during ageing. For example, we showed that this phenomenon is due to crack propagation occurring at the beginning of ageing. A finite element model has then been developed to study more accurately the propagation mode of these cracks. Thermal stresses are likely to propagate cracks through the film in pure opening mode. When cracks are long enough, thermal stresses can generate propagation parallel to the interface. This would result in the decrease of the measured adhesion with a possibility of flaking.

INTRODUCTION

In the space vacuum, thermal exchanges of a satellite with its environment are only radiative. Surface treatments are often used in space applications to provide specific thermo-optical characteristics to materials [START_REF] Shrestha | Microstructural and thermo-optical properties of black Keronite PEO coating on aluminium alloy AA7075 for spacecraft materials applications[END_REF][START_REF] Mccroskey | Aeroglaze Z306 black paint for cryogenic telescope use: outgassing and water vapor regain[END_REF] and thus, to passively manage thermal control. Black inorganic anodized aluminium alloys are then used for their high solar absorptivity and normal emissivity (α s >0.93; ε n >0.9), [START_REF] Goueffon | Black anodic coatings for space applications: Study of the process parameters, characteristics and mechanical properties[END_REF][START_REF] Sharma | Anodizing and inorganic black coloring of aluminum alloys for space applications[END_REF][START_REF] Levesque | Black anodize as a thermal control coating for space station freedom[END_REF]. In addition, those black coatings are especially suitable near optical instruments to avoid stray lights. Moreover, the use of inorganic dyes limits risks of contaminating the spacecraft by outgassing.

With the successive passages of a satellite in the shadow of the Earth in orbit, a satellite is constantly exposed to thermal cycling under vacuum. The temperatures can typically vary between -140°C and 140°C. Many cases of flaking of black anodic films have been observed after thermal cycles carried out to simulate the space environment [START_REF] Goueffon | Black anodic coatings for space applications: Study of the process parameters, characteristics and mechanical properties[END_REF]. Only anodic films on aluminium alloys from series 2XXX and 7XXX are affected by this anomaly. Actually, in orbit, such particles could contaminate instruments of the satellite and reduce its lifetime.

The black anodizing process and the influence of its parameters on the film morphology were studied in a previous paper [START_REF] Goueffon | Chemical characteristics, mechanical and thermo-optical properties of black anodic films prepared on 7175 aluminium alloy for space applications[END_REF]. It has especially been shown that the Young modulus of black anodic films is directly depending on the initial porosity and thus, on the anodizing temperature. A temperature range varying from 5±0.5°C to 25±0.5°C is actually leading to increasing porosities from 10±2% to 50±4% and a decreasing Young modulus from 77±3GPa to less than 10GPa. In this context, the aim of this paper is to understand the mechanisms leading to flaking occurring during thermal cycling. Two types of black anodic films were realised with anodizing temperatures of 20°C or 25°C. Indeed, cases of flaking were observed on samples anodised at ambient temperature. Those samples were first submitted to thermal cycling and the evolution of their adhesion on the aluminium alloy substrate was then evaluated by two methods: peel and scratch-tests. The influence of the cycling conditions (temperatures, number of cycles) on films behaviour was particularly studied. A finite element model was developed to more accurately understand the mechanisms of crack propagation during thermal cycling leading to flaking. The influence of the hot temperatures (>20°C) and cold ones (<20°C) on the stress state in the film was particularly investigated.

EXPERIMENTAL PROCEDURE

The substrate material was the 7175T7351 aluminium alloy, often used in space industry. Its chemical composition in weight percent is: 1.62%Cu, 0.12%Fe, 2.42%Mg, 0.01%Mn, 0.06%Si, 5.75%Zn, 0.041%Ti, 0.21%Cr and Al the remainder.

Black anodizing process

The process of black anodizing used followed the ESA Standard [START_REF]ESA ECSS-Q-70-03A, Black anodizing of metals with inorganic dyes[END_REF] for spacecraft design involving four main consecutive steps: surface pretreatments, anodizing, coloring and sealing. Samples were rinsed with distilled water after each step. The alloy sheet (3x20x40mm) was degreased with ethanol and then etched in Na 2 CO 3 and Na 3 PO 4 aqueous mixed solution for 5min at 93°C and neutralized in aqueous HNO 3 for 3min at room temperature. The aluminium sheet was used as anode and a lead plate as counter-electrode in the electrochemical cell. The anodizing was run for 60min in the galvanostatic mode (Ja=1.2±0.1A/dm²) using a sulfuric acid solution (150g/L) thermally regulated. The electrolyte temperature is not clearly defined in the specification and is generally the room one. Two different electrolyte temperatures were then used: 20°C and 25°C resulting in oxide films with initial porosities of 40% and 50% respectively [START_REF] Goueffon | Black anodic coatings for space applications: Study of the process parameters, characteristics and mechanical properties[END_REF]. The anodized part was colored by precipitation of black cobalt sulphide (CoS) in the pores [START_REF] Scherer | Optimization and microstructural analysis of cobalt-sulphide-pigmented aluminum oxide[END_REF]. It consists in an immersion in two successive baths: firstly for 15min. in a solution of cobalt acetate (200g/L) regulated at 43±2°C, and then for 10min in a bath of ammonium hydrosulphide (30g/L) at room temperature. The sealing step occurred for 25min in a bath regulated at 98±2°C and made up of nickel acetate (NiCH3COO, 4H2O) and boric acid (5g/L both).

Thermal cycling

Space environment has been simulated in an environmental chamber Sun Electronic System EC11. Thermal cycles were performed in a nitrogenous atmosphere between -80°C and 80°C. A large number of parts of a satellite are in this temperature range in space applications. The heating and cooling speed are 10°C per minute with dwell times of 5 minutes. Each cycle starts by the heating of the samples at 80°C before cooling.

Adhesion measurement

The peel-test was first used to evaluate adhesion of the films. A 2cm width aluminium tape was stuck on films. A tensile test device was used to remove the tape with a constant speed of 500 mm/s. the removal force is perpendicular to the surface of the film. If any particle is detached from the film, the sample is considered as failed.

The adhesion of the films was also evaluated using a scratch-test device (CSM Revetest instrument) with a diamond stylus (Rockwell, 200µm radius tip). Scratchtests were configured with an increasing normal load (loading speed: 30N/min; advance speed: 5mm/min). The resulting scratch-print of 5mm length was observed by scanning electron microscopy (SEM Philips XL 30 ESEM) to determine the corresponding normal load. The failure criterion was the flaking of the black anodic film along the scratch. A critical load was defined as the minimum normal load applied on the stylus leading to the lateral flaking of the black anodic film. Given values are averages of six different scratches on each sample.

Finite element model

A finite element simulation was used to evaluate mechanical behavior of a thin film on its substrate. SAMCEF was used to perform isotropic thermo-elastic calculations with 2-dimensional plane strain hypotheses. The geometry used was a substrate (thickness of 3 mm) with an anodic film of 20 µm. Different geometries of cracks in the film have been used and will be described below. Mesh size has been validated by comparing K I with theoretical values in the simple case of a straight crack perpendicular to the surface.

Materials

The mechanical properties of the substrate are well known (E s = 72 GPa, ν s = 0.33, α s = 23.4 10 -6 K -1 ). The properties of black anodic films were studied in a previous paper [START_REF] Goueffon | Chemical characteristics, mechanical and thermo-optical properties of black anodic films prepared on 7175 aluminium alloy for space applications[END_REF] and are assumed to be (ν f = 0.28, α f = 13 10 -6 K -1 ). The young modulus range simulated is 10 GPa / 77 GPa and covers a porosity range of 10% / 50%. The young modulus of those two types of samples experimentally aged was difficult to accurately evaluate by nanoindentation but seems a bit lower than 10 GPa [START_REF] Goueffon | Chemical characteristics, mechanical and thermo-optical properties of black anodic films prepared on 7175 aluminium alloy for space applications[END_REF].

Loadings

The loading is only thermal, constant and uniform. The temperature chosen are either -80°C or 80°C while the reference temperature of the system is 20°C. All stresses in the model are then due to differential thermal dilatations between the film and the substrate.

RESULTS AND DISCUSSION

Characterisation of the films before ageing

The thickness of the films is approximately 20µm. SEM observations of the black films presented on Fig. 1 show the presence of microcracks. The difference of coefficients of thermal expansion between the film and the substrate lead to the crazing of the film during the sealing step at 98±2°C [START_REF] Goueffon | Mechanical behavior of black anodic films on 7175 aluminium alloy for space applications[END_REF][START_REF] Zhao | The cracking behaviours of anodic films on 1050 and 2024 aluminum alloys after heating up to 300°C[END_REF]. In addition Liu et al. [START_REF] Liu | The effects of sealing on cracking tendency of anodic films[END_REF] showed on films isolated from their substrates that the precipitation occurring during the sealing step can induce internal stresses leading to the formation of cracks.

Fig. 1. SEM surface observations of two black anodic films with initial porosities of 40% and 50%

The cracks are wider, and then probably deeper [START_REF] Zhao | The cracking behaviours of anodic films on 1050 and 2024 aluminum alloys after heating up to 300°C[END_REF], on the samples with 50% of initial porosity. This is explained by the decrease of the limit tensile stress of films when the porosity is increasing especially for porosities higher than 30% [START_REF] Goueffon | Mechanical behavior of black anodic films on 7175 aluminium alloy for space applications[END_REF].

Before ageing, in spite of this crazing, peel-tests were successful on both types of samples: no particle was detached by the tape. There is no anomaly detected at this step. The normal load measured by scratch-test is then about 20 N and 17 N for the initial porosities of 40% and 50% respectively.

Adhesion after thermal cycling

The adhesion of those two types of black anodic films was measured after cycles between -80°C and 80°C.

Tab.1. Results of peel-tests after cycles between -80°C and 80°C successful (/) or failed (X) The results of peel-tests performed after ageing are presented in Tab1. Whatever the number of cycles, particles were detached by the tape from samples with 50% initial porosity (Fig. 2). Tests performed on the other type of anodic films are successful except after 3 and 50 cycles where several particles (< 1 cm -2 ) where detached.

Fig. 2. Surface SEM observation of a black anodic film after a failed peel-test

These observations are coherent with scratch-tests results. The adhesion of films with 40% of initial porosity is not affected by this ageing, whatever the number of cycles performed (Fig. 3). Nevertheless, a decrease of the adhesion is detected on the other type of samples as soon as the first cycles.

Fig.3. Critical load evolution with the number of thermal cycling performed between -80°C and 80°C

The decrease of adhesion during thermal cycling is associated to crack propagation under stresses. Two types of stresses can be identified:

-First, the coefficient of thermal expansion of anodic films has been evaluated to 13.10 -6 K -1 [12, 13] while the aluminium alloy's one is about 23.4.10 -6 K -1 . This difference of thermal expansion results in tensile stresses in the film when T > 20°C and compressive ones when T < 20°C.

-Secondly, a dehydration of the film may occur during thermal cycling. It may partially relax the internal compressive stresses of the film due to hydration during the sealing step [14].

These stresses may be able to cause the growth of cracks, reducing the load to apply to pull off particles. The difference of behaviour between the two films can be attributed to their difference of limit tensile stress. The cracking resistance is lower in the film with the highest initial porosity [START_REF] Goueffon | Mechanical behavior of black anodic films on 7175 aluminium alloy for space applications[END_REF]. Actually, there is crack propagation if its stress intensity factor (K) becomes 40% 50% 100 µm higher than the critical stress intensity factor (K c ). K c is only depending on the materials but difficult to evaluate in the case of anodic films.

Model of crack growth through the film

The SEM observations of the fractures induced by peel-test or scratch-test coupled with EDX analysis show that the fractures are cohesive in the film. Actually, the anodizing is an electrochemical conversion of the alloy surface. The interface between the film and the substrate has always been considered strong [15]. A finite element simulation has been used to understand the crack growth mechanisms within the film.

Geometry

A 2-D model of the cross-section of a thin film (20 µm) on its substrate (3 mm) before ageing has been created. The influence of a straight crack perpendicular to the surface with a width of 1 µm and a length varying from 9 to 18 µm on the stress state was observed (Fig4).

Model exploitation

Heating the sample at 80°C generates tensile thermal stresses in the film depending of its young modulus. The crack is then solicited in pure opening mode (Mode I). The evolution of K I with the young modulus is given on Fig. 5 for different crack lengths.

The strains of the film are imposed by the substrate. Thus, an increase of the film's Young modulus results in higher stresses. When the young modulus is close to the substrate's one (72 GPa), the stress intensity factor (SIF) is increasing with the length of the crack as it is the case in an homogenous material.

When the modulus of the film is much lower, the stress intensity factor becomes independent of the crack length. This can be associated to a screen effect of the substrate. When the crack approaches the interface, strains are blocked by the more rigid substrate. That limits the increase of the SIF with the length of the crack.

Fig. 5. Evolution of the stress intensity of the crack at 80°C with the young modulus of the film for different crack lengths Whatever the modulus of the film, cracks perpendicular to the surface can grow under thermal stresses for temperatures higher than 20°C.

Lower temperatures than 20°C result in a stress intensity factor equal to zero; compressive stresses close the crack.

Model of crack growth parallel to the interface

A crack, parallel to the interface may propagate to generate flaking. A numerical model was built to study the mechanisms of growth of these cracks. They can be initiated by the presence of bubbles or defects in the anodic films.

Geometry

A kinked crack was modelled with a varying depth as presented on Fig6. The width of the crack is 1µm and we have shown that the radius of curvature in the kinked part of the crack has no influence on the stress state near the crack tip, in the area used to calculate the stress intensity factor. Film Substrate Fig. 6. Model (on the left) of the geometry of a kinked crack in a film compared with a SEM cross-section observation of a film with 40% of initial porosity (on the right)

Model exploitation

The differential thermal expansions are generating tensile stresses in anodic films at 80°C. The presence of a kinked part is locally modifying the stress state. Tensile stresses are present on one side of the crack while compressive ones are on the other side (Fig. 7). The crack is then solicited by the sliding mode (Mode II). The stress intensity factor associated to this mode (K II ) was calculated (Fig. 8).

Fig. 8. Evolution of K I and K II with the young modulus of the film at -80°C and 80°C for crack depth of 17 µm At 80°C, the crack is solicited in a mixed mode (I and II) while there is a pure mode II at -80°C. The value of K II is mainly depending on the temperature difference with the reference (20°C), thus K II (-40°C) ≈ K II (80°C). That explains that the value of K II at -80°C is much higher than at 80°C. Both hot and cold parts of thermal cycles are thus able to lead to crack propagation parallel to the interface.

Fig. 9 shows that the maximum stress intensity (K I and K II ) factor is obtained when the crack kink is the deeper in the film, whatever the considered temperature.

CONCLUSION

Before any ageing, no adhesion anomaly is detected on black anodic films, in spite of the presence of cracks due to the process. Adhesion of these films can be affected by thermal cycling. When it is the case, peel-test and scratch-test reveal the anomaly after the first cycle.

The anodizing temperature seems to be a key parameter for flaking problems. Actually, an increase of 5 degrees changes the initial porosity of the film from 40% to 50%. As a result, black films anodized at 20°C are not visibly affected by ageing between -80°C and 80°C while films realised at 25°C are flaking.

Numerical simulation has shown that hot temperatures can generate crack propagation through the thickness of the film in opening mode. If the crack is deep enough, a defect in the film (oxygen bubbles for example) can initiate a cracking near the interface film/substrate. The propagation along the interface could then occur under pure mode II or under mixed mode. It would result in the detachment of particles and possibly in the pollution of particles.

Before ageing, cracks in films with 50% of initial porosity are wider and probably deeper than in the film with 40% of initial porosity. The small difference of young modulus between those two types of films is not explaining the difference of behaviour observed during ageing. The more porous films may have a lower critical stress intensity factor leading to easy crack propagation. Thus, thermal stresses between -80°C and +80°C may then be sufficient to cause propagation along the interface only in the films with the highest porosity, resulting in the decrease of adhesion observed.

Fig. 4 .

 4 Fig.4. Geometry of the crack in the numerical model (on the right, visualisation of σ xx ) compared with a crack observed by SEM on a black anodic film with 40% of initial porosity before ageing (on the left)

Fig. 7 .

 7 Fig.7. Visualisation of σ xx in a film with a kinked crack at 80°C

Fig. 9 .SubstrateFig. 10 .

 910 Fig.9. Evolution of the stress intensity factor in a film (E f = 10 GPa) with the kink's crack depth.Thus, defects (bubbles) near the interface are more likely to propagate under thermal stresses. Hot temperature (> 20°C) would lead to propagation under mixed mode while cold temperatures (< 20°C) are loading the crack with pure mode II. Crack propagation could lead to a decrease of the critical load measured by scratch-test or the failure of the peel-test. Fig.10shows a cross section of a black anodic film (50%

AKNOWLEDGEMENT

The authors thank Gregory Aldebert for his participation to this work.

on 2024 aluminum alloy after heating up to 300 °C, Surf. Coat. Tech., 203, 9, 2009, 1244-1251. 12. Zhou J., et al.