Sampling of A. thaliana in Pyrenees mountains allowed to identify 30 independent locations at different altitudinal levels all along the Pyrenees mountains (Table 1). These 30 independent sampled locations have been mentioned as natural populations and corresponded to 341 individuals. Their taxonomic belonging to A. thaliana species was confirmed through DNA sequencing (cf. genetic analysis). Population names correspond to the first four letters of the closest village or location where the plants were found (Table 1). One locale population Lant (Lanta, 250 m, France, Bartoli et al., 2018), Col (Columbia, 200 m; Poland) and Sha (Shahdara, 3,400 m; Tajikistan) were used as out-groups. Col and Sha have already been studied in the lab under the same growth conditions (Duruflé et al., 2017). Field-collected seeds from about 10 individuals per sampled locations (Table 2) were reproduced one time in a growth chamber to obtain homogeneous batches of seeds and to reduce maternal effects, before phenotyping. DNA was extracted from all individuals, and only one individual was selected per sampled locations for phenotyping. Seeds were sown in Jiffy-7^® peat pellets (Jiffy International, Kristiansand, Norway). After 48 h of stratification at 4°C in darkness, plants were grown at 90 μmol photons m^-2 s^-1 light intensity, 70% humidity and 16 h light/8 h dark photoperiod at two different temperatures: 22°C (optimal growth condition for Col) or 15°C (sub-optimal temperature for Col but potentially optimal for some of the Pyrenean populations). For rosette phenotyping and biochemical analyses, four- and six week-old rosettes were collected, respectively, at 22 and 15°C, corresponding to emergence of the first flower buds for Col and Sha (Boyes et al., 2001). Each phenotyping was performed in triplicate with 2 to 5 plants per replicate. Plants were grouped by population but their positions were randomized between the three replicates and during the experiments.

Climatic variables were obtained from WorldClim dataset¹ (Hijmans et al., 2005). The values used are the mean of 30 years (1960 to 1990) with a resolution of ca. 1 km² per grid cell. Solar radiation reading were obtained from the Photovoltaic Geographical Information System of the European Communities (Huld et al., 2012) for the 2001–2012 period². The specific solar radiation of Sha was estimated by linear regression of altitude level on solar radiation of the Pyrenees mountains.

Genomic DNA was extracted using a standard CTAB protocol from leaves. Five pairs of primers (sequences given in Supplementary Figure S1A) were designed to amplify 5 genomic areas of about 1000 bp each. The 5 regions are distributed on the 5 chromosomes of A. thaliana; they encompass ca. 500 bp upstream of the ATG and ca. 500 bp both coding and non-coding region (including exon and intron). They correspond to the loci of 5 CIII Prx [At1g44970 (AtPrx09), At2g41480 (AtPrx25), At3g50990 (AtPrx36), At4g33870 (AtPrx48), At5g39580 (AtPrx62)]. PCR were performed using high fidelity recombinant Pfu DNA polymerase (Promega, Madison, WI, United States) according to the manufacturer’s instructions. PCR products were sequenced using the same primers and on both strands to ensure reading quality. Paired end sequencing and manual editing using BioEdit (Hall, 2011) allowed a good sequencing quality. The five sequences obtained were concatenated and used for genetic and phylogenetic analysis. In order to confirm that a subset of populations was homogeneous (Camu, Gava, Grip, Hern, Hosp, Jaco, Pont, Prad and Roch), PCR amplification of 3′ETS 45S rRNA genes was done using primers flanking the 3′ETS variable region. The size of the product depends on the number of repeats and corresponds to variant of 45S rRNA genes (Pontvianne et al., 2010). The QGIS software 2.18³ was used to georeference the A. thaliana populations on a map and to display associated pie charts illustrating the proportion of 45S rDNA genotypes identified for each population.

The sequences of 30 Pyrenean populations, the 3 out-group (Lant, Col and Sha) and 22 accessions available from the 1001 Genomes Project (11 in South France and 11 in North Spain the most closely to the Pyrenees mountains) were analyzed (Supplementary Tables S1, S2). Tajima’s D neutrality test was performed using DnaSP v5 software (Librado and Rozas, 2009). Tajima’s D was calculated using the concatenated sequence of the 5 CIII Prx loci. It allows to compare the total number of segregating sites (polymorphisms) to the average pair-wise differences between sequences (Tajima, 1989). This test distinguishes DNA sequences evolving under selective pressures from those evolving neutrally. For Col, Sha and the 22 French/Spanish accessions from the 1001 Genomes Project, only one individual was used to minimize the impact of the frequencies of these new haplotypes but also because no information regarding the genetic diversity of these accessions was available. Across the 5 loci, the haplotype class of each individual sequence (i.e., the individuals were homozygotes) has been determined using the Col sequence as reference. Briefly, each nucleotide polymorphism was encoded as binary (i.e., 0/1, with 0 being the “Col” allele) and according to this matrix, the haplotypes have been defined when they possessed at least one polymorphism relative to Col. Population structure was analyzed using a Bayesian clustering method implemented in the STRUCTURE software version 2.3.4. (Pritchard et al., 2000). Assuming K differentiated genetic clusters in the sample; this method allows estimating the proportion of membership of any individual to any of the K clusters. For each run with a K value, 10⁶ iterations and a 10⁵ burn-in period options were chosen. For each value of K (between 1 and 10), 10 independent runs were performed, and likelihood values obtained from these 10 runs were averaged. To identify the appropriate K value, Evanno’s method based on data likelihood variation over successive K was used (Evanno et al., 2005). The optimal K value was estimated as the one with the largest delta K value (Supplementary Figures S2A,B).

The analysis involved 375 nucleotide sequences (from 341 natural individuals from Pyrenean populations, 10 individuals from Lant, the 2 references Col and Sha and 22 accessions from the 1001 Genomes Project). Phylogenies can be inferred separately from each gene or directly from concatenated supergene with various advantages (Gadagkar et al., 2005; Kubatko and Degnan, 2007). The concatenation approach has been chosen, and the 5 CIII Prx sequences were assembled for a total of 4074 positions in the final dataset. After alignment, all positions containing gaps and missing data were eliminated. A phylogenetic tree was inferred with concatenated sequences by using the maximum likelihood method based on the Tamura-Nei model (Tamura and Nei, 1993). 100 bootstraps were generated to evaluate phylogenetic robustness. These analyses were conducted with MEGA v6 (Tamura et al., 2013).

Most data analyses were performed with R software (version 3.2.3). Student’s statistical tests for two samples were carried out in order to determine the temperature effects on the populations. In addition, Duncan’s multiple range tests were performed with the package agricolae⁴ to allow a better visualization of the different sets when they correspond to averages. Kruskal–Wallis and Wilcoxon non-parametric tests were carried out to compare groups. To investigate the underlying variation between population data (phenotypic), Principal Component Analysis (PCA) with multilevel function was used with the mixOmics package (Lê Cao et al., 2011) available in CRAN⁵. This combination of multilevel and multivariate methods allowed distinguishing variations between populations and between temperatures. Multilevel PCA was done only with contrasted data collected at 15 and 22°C. In addition, distance between the same populations according to the temperature in the scaled PCA multilevel was measured. Graphical representation of the STRUCTURE results was done with the “leaflet” package⁶.

The 30 Pyrenean populations of A. thaliana as well as Lant, Sha and Col, were described using the altitude (m) and the following 5 climatic variables (Table 1): annual minimum, annual maximum and annual mean temperatures (°C), total annual precipitation (mm) and total annual UV radiations (kWh/m²).

The variables describing temperature spectrum (minimum, mean and maximum) are highly positively correlated but negatively correlated with altitude, annual rain accumulation and solar radiation (Supplementary Figure S4A). To analyze the environmental datasets (climates and altitudes), PCA were used to calculate principal components (PCs) that explain the climatic variance (Supplementary Figure S4). The PCs score approach is routinely used for climate quantification (Wolfe and Tonsor, 2014) and represent an index that well describes natural environmental conditions of each population (Table 1). In the PCA, the climate PC1 explained 74% of the multivariate variance across the six variables, while climate PC2 only explained 14% (Supplementary Figure S4B). Climate PC1 is strongly associated with most variables, while climate PC2 is only associated with the atypical climatic profile of Sha (Supplementary Figure S4A). Thus, the environmental data characteristics of each population can be summarized with climate PC1 value (Table 1) and will be used for further integrative analysis.

In this study, the polymorphism of five CIII Prx loci (AtPrx09, AtPrx25, AtPrx36, AtPrx48 and AtPrx62) was analyzed. Those loci played major roles in Col, during most of development processes and in response to many biotic and abiotic stresses (Francoz et al., 2015) without a priori implication in altitude adaptation, making them good markers for population genetic structure analysis of each geographical populations. These five amplified sequences include both coding and non-coding regions (Supplementary Figure S1B). They have been selected for their significant degree of polymorphism [Single Nucleotide Polymorphism (SNP) and Insertion/deletion (Indel)] detected between accession sequences from 1001 Genomes Project (Weigel and Mott, 2009).

In addition to the sequences obtained from the 341 Pyrenean individuals, the sequences from populations used as out-groups (10 individuals of Lant, and one Col and Sha) and those of 22 accessions available from the 1001 Genomes Project (1001 Genomes Consortium, 2016) and localized in the south of France and the north of Spain were added. These accessions originating from both sides of the Pyrenees mountains improved the sampling in a regional context. Based on the Col sequence, 338 positions of SNP and Indel were identified and used to define different haplotypes (combination of alleles of different markers in the sequences). Tajima’s D was calculated and it was not significantly different from zero (D = 0.420, p > 0.10), meaning that null hypothesis of neutral mutation drift equilibrium for these sequences could not be rejected. Consequently, these markers could then be considered as neutral and used for genetic analysis of population structure. The polymorphism within the populations was evaluated using different parameters: the haplotype diversity and frequency distribution, the number of SNPs and the percentage of identity of haplotypes contents in each population (Supplementary Figure S5, Supplementary Data Sheet S1, and Table 2). This analysis allows to determine that 9 of the 30 Pyrenean populations (i.e., Camu, Gava, Grip, Hern, Hosp, Jaco, Pont, Prad and Roch) are homogeneous (i.e., with few intra-population variation, see Table 2).

To analyze the genetic structure and polymorphism of the Pyrenean A. thaliana populations, two different approaches were used: a Bayesian clustering analysis and a phylogenetic analysis. These two analyses complement to evaluate the different polymorphism parameters in the genetic analysis.

In order to infer population structure, haplotype frequencies were used to analyze the genetic diversity within and between populations. In the clustering analysis implemented in STRUCTURE, it is possible to test various genetic clusters (K value), each characterized by a set of haplotype frequencies. These STRUCTURE analyses allowed to identify two genetic clusters for the Pyrenean populations (Supplementary Figure S2). Populations with a membership proportion of >0.8 to a cluster are considered as homogeneous (non-admixed) and populations with <0.8 proportion are considered as heterogeneous (admixed) (Fukunaga et al., 2005; Uzun et al., 2014). The clustering analysis identified 36 populations considered homogeneous among the 55 analyzed. In this Bayesian clustering analysis, the accessions from 1001 Genomes, the populations used as out-group (Lant, Col and Sha) and 63% of the Pyrenees mountains populations were found in the same genetic cluster (Supplementary Figure S2C). From all the Pyrenees populations, 11 (37%) shared the same genetic cluster different from the known French and Spanish populations (Figure 1). Three genetic clusters were defined with the Bayesian clustering analysis as the following: (i) populations assigned to clusters 1 (green, Pyrenean-specific) and 2 (pink) are homogeneous populations with membership proportion > 0.8, (ii) populations assigned to cluster 3 are considered as population of mixed ancestry because they do not belong to any of the two previous clusters.

A phylogenetic analysis was performed to obtain a complementary view of the natural genetic diversity among the Pyrenean A. thaliana populations. Genetic relationships among the different haplotypes were determined by the Neighbor-Joining (NJ) method using the concatenated sequences of the 5 CIII Prx markers. The tree segregated the populations into three major clusters (Figure 2). Moreover, several populations were grouped on the same subtree (e.g., Hosp, Hern) or close to each other on the same branch (e.g., Roch, Grip, Jaco) as expected based on polymorphism parameters (Table 2).

The germination event is crucial for the adaptation of plant populations to different climatic conditions. To test seeds viability of individuals of the 30 Pyrenean populations, their germination rate was quantified at 22 and 15°C. Only one individual per population was randomly chosen regardless the results of the genetic analyses. It should be noted that the Mong population was not included in phenotyping as its natural location was destroyed by human activity and seeds were not hence collected.

Molecular phylogenetic analysis (maximum likelihood method) based on five independent peroxidase gene markers of 341 individuals from the 30 A. thaliana Pyrenean populations, 22 accessions taken from the 1001 Genomes Project (11 from the South of France and 11 from the North of Spain) and the 3 out-groups Lant, Col and Sha (10 individuals for each). Color represents the population allocations into their estimate membership proportions in each genetics cluster determined by STRUCTURE and illustrated in the Figure 1. Numbers in brackets stand for the number of individual per cluster (collapsed branches) and numbers without brackets represent the individual within the population. Percentage of bootstrap values higher than 50 are indicated above the main branches. The branches and clusters are colored with the same colors as STRUCTURE analysis (Figure 1) with the admixed populations in gray.

The germination rate was constant for all the Pyrenean populations tested with a rate higher than 90% (data not shown). It demonstrated a high seeds viability in all the populations independently of the temperatures used for the germination assays. However, the germination speed differed between the populations.

In A. thaliana, successive early steps of germination comprise testa rupture (TR) and endosperm rupture (ER) (Supplementary Figure S3a). The percentages of TR and ER were evaluated at 15 and 22°C, 24, 30, and 48 h after seed imbibition (Figure 4 and Supplementary Figure S3b). The germination speed was largely reduced at 15°C and varied considerably between populations. For example, Hosp, Col, Jaco, Mere, Arag and Guch had high percentages of TR and ER at 22 and 15°C. Whereas Pont, Gava, Savi, Grip, Prat, and Fos germinated significantly slowly at 22 and 15°C (Supplementary Figure S3a). Interestingly, it should be noted that 4 of the 6 most slowly germinating populations at both temperatures, originate from high altitudes. The TR/ER ratio was significantly higher for Eaux, Herr, Biel, Gedr, Bedo as compared to the other Pyrenean populations, meaning that in the former populations the testa layer could be weaker or that the embryo growth may be slower than in the majority of the other accessions. Conversely, Hosp, Arag, Mere, Jaco and Guch showed a profile similar to that of Col, with a percentage of ER higher than the TR (Figure 4). This could be explained by the resistance of the testa or the different forces (e.g., pressure, turgescence) applied by the embryo (Müller et al., 2006).

The individuals, used for germination assays, were also phenotyped at later developmental stages. To better understand phenotypic variations in response to temperature variation, several traits of A. thaliana grown at 22°C (optimal growth condition for Col) were compared to those of plants grown at 15°C (to reproduce low temperature of altitudinal growth conditions). Contrasted rosettes phenotypes were observed between the two temperatures (Figure 5A,B) and quantified such as the number of leaves (Supplementary Figure S6), the rosette weight (Supplementary Figure S7) and their diameter (Supplementary Figure S8), the time to reach the bolting stage (Supplementary Figure S9), the rosette anthocyanin (Supplementary Figure S10) and chlorophyll contents (Supplementary Figure S11). In parallel, the floral stem phenotypes were quantified such as the number of cauline leaves on the main stem (Supplementary Figure S12), the stem length (Supplementary Figure S13) and diameter (Supplementary Figure S14), and the stem growth speed (Supplementary Figure S15). It should be noted that Pont appeared to be a “winter-annual” Arabidopsis population (Chouard, 1960; Michaels and Amasino, 1999; Gazzani et al., 2003). Then only the rosette data were recorded for Pont because floral stems data have been obtained after transfer to 4°C for 4 weeks, they were too different to be compared with those of the other Pyrenean populations.

In general, the number of leaves, the weight, and the diameter of the rosettes were higher at 15°C than at 22°C, except for Bran and Col and Sha whose rosettes had a smaller diameter at 15°C, (Duruflé et al., 2017) or Pont which had a reduced rosette weight at 15°C (Supplementary Figures S6A, 7A).

In parallel, increased stem diameter, stems length, and cauline leaves number were observed in plants grown at 15°C, with the notable exceptions of Hern and Prad (less cauline leaves), and Biel (shorter stems at 22°C than at 15°C). The stems growth speed was also significantly faster at 15°C. For example, Eaux, Fos and Chau stems grew slower at 22°C than at 15°C (Supplementary Figure S15).

To evaluate the level of resistance to stress of A. thaliana Pyrenean plants grown at different temperatures, anthocyanin content was quantified in the leaves of individuals cultivated at 15 and 22°C. Most of the populations (e.g., Argu, Jaco, Eaux and Bier) showed an increase of anthocyanin content of rosettes grown at 15°C as compared to 22°C (Supplementary Figure S10).

The Pyrenean populations could be divided into two groups based on their relative chlorophyll content when grown at 15 or 22°C: those with significantly more chlorophyll at 15°C than at 22°C (e.g., Mere, Chau, Eget, Bran, Biel) and those with more chlorophyll at 22°C than at 15°C (e.g., Cast, Hosp, Prad, Roch, Hern). It should be noted that regardless of the variation between the two growth temperatures, all the populations revealed about the same chlorophyll content at 15°C (Supplementary Figure S11D).

All the phenotypes could have been analyzed independently to study the acclimation mechanisms of the Pyrenean populations of Arabidopsis to the different growth temperatures. In fact, the 10 traits measured (e.g., number of leaves, anthocyanin content) varied not only between the growth temperatures but also between the populations/accessions; hence, PCA was performed to get more integrative results. Multilevel methods of PCA can be used if repeated measurements from different treatments are applied on the same individual. This multilevel approach was developed to highlight treatment effects within subject separately from the biological variation between subjects (Liquet et al., 2012). It was dedicated to take into account potential artifact in the data due to repeated measurements and allows investigation of the between subject variation which is completely separated from the within subject variation (Westerhuis et al., 2010). This new approach, not widely used in the phenotypic studies, is very efficient to study the effect of different conditions within subjects separately from the variation between populations. Since our dataset fulfilled these conditions, a multilevel PCA of the phenotypic data was performed (Figure 6A,B).

According to this multilevel analysis, populations were clearly separated by the temperature effect in the first axis (55%). This separation was mainly due to the chlorophyll overall content which was more abundant in rosettes from plants cultivated at 22°C as compared to 15°C, and to the other phenotypic parameters generally higher at 15°C than at 22°C in most of the populations. Indeed, highland population such as Prad and Hosp appeared to contain more chlorophyll at 22°C than at 15°C unlike to the lowland populations Chau and Biel. The second axis differentiated Pyrenean populations according to their speed to reach the bolting stage, which was negatively correlated with the floral stem speed regardless of the temperature. For example, Bran, Fos and Eaux reach the bolting stage earlier when grown at 22 vs. 15°C and have a faster stem growth speed than the other Pyrenean populations when cultivated at 22°C. Therefore, the time to reach bolting stage and stem growth speed appear to be independent and not correlated with the developmental stages, and hence useful to compare different genotypes.

The effect of genetic clusters or environmental data on phenotypic trait variation was evaluated by Kruskal–Wallis non-parametric tests Figure 7A,B). Only few traits such as chlorophyll content, stem diameter are correlated with genetic and the environmental groups.

In the populations assigned to the three distinct genetic clusters by STRUCTURE analysis, rosettes and floral stems phenotypes were impacted by growth temperature. Indeed, 9 of the Pyrenean populations considered as homogeneous (cluster 1 and 2) have systematically wider stems than the population of mixed ancestry (cluster 3) at both growth temperature conditions (Figure 7C). Moreover, rosettes of populations from the cluster 3 are smaller than those of populations from cluster 1 at the same temperature. At 22°C, the same tendencies are visible, although not significant despite the low variability of the data in this genetic cluster (Figure 7D). Because populations of cluster 1 are specific of the Pyrenean mountain, one phenotypic characteristics of A. thaliana in the Pyrenees is to have bigger rosette at low temperature growth conditions.

Finally, climate PC1 value obtained from normed PCA of six environmental datasets, allows classifying populations as a function of their environmental conditions considering altitude levels. The populations were distributed as plain (PC1 value > 1), medium (0 < PC1 value < 1) and high altitudes (PC1 value < 0) (Table 1 and Figure 7A). Some correlations exist between phenotypic traits and the original environmental conditions characterized by the climate PC1 value. For example, populations living at high altitude showed significantly higher chlorophyll content at 15°C compared to the other populations. However, this was not observed at 22°C (Figure 7E).

All the climatic variables are usually well correlated with altitude, although the field topology may play an important role. Precipitation and seasonality exert the largest influence on the regional climate variation (Körner, 2007). Pyrenees mountains have a temperate climate where the gradient of precipitation is dependent on altitude. It is difficult to compare the characteristics of all populations at once but several conclusions can nevertheless be drawn from environmental data. Climate PC1 value obtained from normed PCA of six environmental datasets, allows classifying populations as a function of their environmental conditions considering altitude levels: population of plain, medium and high altitudes. The climate PC1 values are highly correlated with the population altitude but some differences could be noted. Eaux has a lower climate PC1 value than other populations living at the same altitude such as Savi or Bier. This could be explained by the existence of lower temperatures at this location compared to those where Hern and Roch grow. On the other hand, although being populations from low altitude, Bran and Gedr have a negative climate PC1 value, as do the Pyrenean populations living above 1100 m of altitude, which could be explained by the lowest temperatures and the high precipitations. Thus, climate PC1 values is useful for comparing environmental conditions and phenotypes data in an integrative way.

Sequence polymorphism is currently used to analyze genetic variability and structure of natural populations (Gulsen et al., 2010; Nemli et al., 2014) and to study evolutionary relationships and genotypic diversity on an intra and inter-specific basis (Gulsen et al., 2010; Uzun et al., 2014). According to our analyses, several geographically distant Pyrenean populations share the same haplotype (Supplementary Table S2). This could mean that these populations have the same origin, either natural or related to anthropogenic introduction/displacement (e.g., seed or soil transportation, road construction…).

According to the 5 CIII Prx markers used, most (70%) of the A. thaliana populations collected in the Pyrenees mountains are genetically separated from the 22 accessions originated from the south of France and the north of Spain and described in the 1001 Genomes Project. The specific structure of these populations in this geographic area makes them very interesting to understand regional population structure in A. thaliana. However, the genetic diversity in the Pyrenees mountains is not just defined by a single genetic lineage; rather, fine-scale patterns of population structures were depicted by the phylogenetic tree. These inter- and intra-population variations may be associated with the variability of the natural environment and may demonstrate a specific adaptation of A. thaliana within this topographic complex area that may have acted as geographic barriers. Nevertheless, they could also have resulted from a prehistoric introduction or from a glacial refuge. It is interesting to note the complementarity of the results obtained with the different approaches (polymorphism analyses, Bayesian clustering and phylogeny). Indeed, populations that are homogeneous according to the Bayesian clustering approach (STRUCTURE) are also those with a high percentage of sequence identity and were found on one branch of the tree (i.e., Grip, Hosp, Hern, Prad and Roch). Conversely, the populations considered as heterogeneous were distributed into two main branches of the tree (i.e., Belc, Fos, Gedr).

Haplotype diversity and frequency (Table 2), Bayesian clustering (Figure 1) or Neighbor-Joining analyses (Figure 2) highlight an unexpected natural and local genetic diversity within and between populations. Satisfyingly, these three approaches identified the same populations as homogeneous (Jaco, Grip, Hosp, Hern, Roch and Prad). Therefore, our study shows that the regional inter-population variability can be important even in highly self-fertilizing plants and must not be neglected in omics projects. On the other hand, the weak intra-population variability that we observed is expected in self-fertilizing plants such as Arabidopsis. Although convincing, conclusions obtained with the markers used about the genetic structure of populations may not be representative of the whole genome and should be supported by more sequencing.

The early step in seed life is crucial for plant survival faced to natural climatic variations. The cell wall composition of the testa could be analyzed to confront these results and to link germination rate to the natural growth environments. The control of germination can be a predominant phenotypic trait for plants to be more adapted to environmental constraints with two different strategies: (i) a rapid germination which allows the plant to grow faster in the frame of favorable climatic conditions and (ii) a slow germination speed which allows the seed to optimize the timing to germinate in different environmental constraints. This later strategy could be correlated to seed dimorphism as observed in Aethionma arabicum (Lenser et al., 2016). In addition, the differential germination speed may explain the contrasted time needed to reach the bolting stage (Toorop et al., 2012) and, consequently increases the chance to grow within diverse climatic constraints. However, as no statistically significant correlation has been shown between germination speed and natural habitats of the populations, no general rule can be found.

Anthocyanin accumulation in leaves is considered as a general stress marker in plants (Mori et al., 2009) but in some situations such as for highland species, they can protect the plant against UV radiations (Close and Beadle, 2003). The anthocyanin accumulation observed for most of the populations grown at 15°C can be a marker reflecting difficulties for these populations to have optimal growth in cooler conditions. However, as expected from its natural localization (1456 m), Pont (and the non-Pyrenean altitudinal control Sha) is well adapted to grow at low temperature regarding its low anthocyanin level even when grown at 15°C.

Chlorophyll content is an indicator for the overall availability of photoreaction centers (Ballottari et al., 2007) as well as of the whole metabolism. The variability of chlorophyll content at low temperature can be due to a compensatory effect to face the reduction of the metabolic activity (Strand et al., 1997) and the variability of Rubisco catalysis present in natura (Stitt and Schulze, 1994; Walker et al., 2013). The natural variability observed between 22°C and 15°C with the Pyrenean populations might be useful to better understand acclimation to multi-climatic changes. Chlorophyll content appears to be a good proxy to characterize the acclimation of populations to more contrasted temperature growth conditions.

In this study, some Pyrenean populations show early bolting stage at 22°C compared to Col and Sha such as Eaux, Cast, Fos, Chau or Gava (Supplementary Figure S9). It appeared that populations with a later bolting stage at 22°C showed a flowering time equal or advanced at 15°C (i.e., Prad, Argu and Belc). In both situations, no relationship was observed with natural habitats of the populations Activity of the regulators of the floral repressors FLC (flowering locus C), FCA and FVE (flowering locus CA and VE) or the SVP (short vegetative phase) could explain this natural variation of flowering time with different ambient temperatures (Simpson and Dean, 2002; Blázquez et al., 2003; Lee et al., 2007).

Contrasted temperature growth conditions highlight different phenotypic variation. As the multilevel PCA (Figure 6C) allows a good separation according to the temperature, the distance separating the two growth conditions of one population can be a good indicator of the global acclimation. A longer distance between the two temperatures means that the phenotypes were more contrasted between the two growth conditions, and therefore suggests a higher phenotypic variation. Two medium and two high Pyrenean populations (Roch, Eaux and Grip, Prad, respectively) showed an elevated phenotypic variation and also a low genetic variability (homogeneous populations). In contrast, the populations from plain and medium environment (Mere, Savi, Chau, Guch, Lant and Bedo) exhibited less phenotypic changes and a large genetic variability. Since, only one individual per population has been phenotyped, it is hazardous to draw general conclusions about phenotypic variation for the whole population. This approach could still help to predict the evolution of plastic responses to environmental changes (van Kleunen and Fischer, 2005; Vitasse et al., 2010; Auge et al., 2017). For example, for a given temperature, wider floral stems of the individuals belonging to homogeneous populations might allow the populations to resist to strong winds in order to sustain seeds dissemination.

These phenotypic analyses showed high intraspecific morphology variations that could be explained by the natural diversity of acclimation to contrasted temperature variations (Hamilton et al., 2015). The absence of clear correlation for all the populations between the rosette weight and the number of rosette leaves could be due to different leaf densities. Indeed, leaves and petioles could be thicker due to a higher quantity of cell wall. Plant cell wall is crucial in growth control, for the structural integrity of the global plant shape (Sasidharan et al., 2011) and contributes to plant architecture variation in response to environmental changes (Le Gall et al., 2015; Houston et al., 2016). These morphological variations may be associated with cell wall modifications (Duruflé et al., 2017). Acclimation of the highland populations to low temperature could be due to these phenotypic characteristics.

Finally, no clear correlation could be established between genetic cluster assignment and the overall environmental data (Figure 7B). Only some phenotypes can be linked to the genetic cluster assignment (floral stem diameter and the anthocyanin content at 15°C) and the environmental data (the floral stem speed at 22°C and the chlorophyll content at 15°C) and could therefore be used as bio-markers for this species. Regarding the genetic clustering, the results exposed highlights the case that populations of mixed ancestry are mostly located at medium altitude (Figure 7A).

In summary, the present characterization of 30 populations of A. thaliana using 5 CIII Prx and rDNA markers allowed to assess the population structure and the genetic relationship between these populations. Genetic structure as well as inter- and intra-population variation emphasized the unexpected variability found in this region. Regarding the various genome sequencing projects (1001 Genomes Project): 9 Pyrenean populations appear genetically close and homogeneous with regard to other accessions known in the region; 12 Pyrenean populations appear genetically close, homogeneous and specific of these mountains and 9 Pyrenean populations appear as populations of mixed ancestry.

This study also revealed phenotypic variations in acclimation of A. thaliana across abiotic gradient modeled here with the temperature change. Some of them are correlated with identified genetic clusters or with environmental data. These analyses contribute to enrich knowledge on abiotic stress acclimation in natural plant populations. Preliminary analyses suggest that the phenotypic acclimation may be due to the cell wall modification affected by low temperature (Duruflé et al., 2017). Transcriptomic, proteomics and cell wall polysaccharides analysis could prove helpful in study cell wall plasticity.

Adaptation of plants to changing environments is a complex process in the context of rapid and transitory environmental changes. It is not possible to conclude that the most adapted phenotype to the temperature elevation is always the most plastic phenotype, but phenotypic trait variation could drive the acclimation of A. thaliana populations to heterogeneous abiotic conditions found in mountain environments.