Phenotypic Trait Variation as a Response to Altitude-Related Constraints in Arabidopsis Populations.

Natural variations help in identifying genetic mechanisms of morphologically and developmentally complex traits. Mountainous habitats provide an altitudinal gradient where one species encounters different abiotic conditions. We report the study of 341 individuals of Arabidopsis thaliana derived from 30 natural populations not belonging to the 1001 genomes, collected at increasing altitudes, between 200 and 1800 m in the Pyrenees. Class III peroxidases and ribosomal RNA sequences were used as markers to determine the putative genetic relationships among these populations along their altitudinal gradient. Using Bayesian-based statistics and phylogenetic analyses, these Pyrenean populations appear with significant divergence from the other regional accessions from 1001 genome (i.e., from north Spain or south France). Individuals of these populations exhibited varying phenotypic changes, when grown at sub-optimal temperature (22 vs. 15°C). These phenotypic variations under controlled conditions reflected intraspecific morphological variations. This study could bring new information regarding the west European population structure of A. thaliana and its phenotypic variations at different temperatures. The integrative analysis combining genetic, phenotypic variation and environmental datasets is used to analyze the acclimation of population in response to temperature changes. Regarding their geographical proximity and environmental diversity, these populations represent a tool of choice for studying plant response to temperature variation.


HIGHLIGHTS
-Studying the natural diversity of A. thaliana in the Pyrenees mountains helps to understand European population structure and to evaluate the phenotypic trait variation in response to climate change.

Natural variations help in identifying genetic mechanisms of morphologically and developmentally complex traits. Mountainous habitats provide an altitudinal gradient where one species encounters different abiotic conditions. We report the study of 341 individuals of Arabidopsis thaliana derived from 30 natural populations not belonging to the 1001 genomes, collected at increasing altitudes, between 200 and 1800 m in the Pyrenees. Class III peroxidases and ribosomal RNA sequences were used as markers to determine the putative genetic relationships among these populations along their altitudinal gradient. Using Bayesian-based statistics and phylogenetic analyses, these Pyrenean populations appear with significant divergence from the other regional accessions from 1001 genome (i.e., from north Spain or south France). Individuals of these populations exhibited varying phenotypic changes, when grown at sub-optimal temperature (22 vs. 15 • C). These phenotypic variations under controlled conditions reflected intraspecific morphological variations. This study could bring new information regarding the west European population structure of A. thaliana and its phenotypic variations at different temperatures. The integrative analysis combining genetic, phenotypic variation and environmental datasets is used to analyze the acclimation of population in response to temperature changes. Regarding their geographical proximity and environmental diversity, these populations represent a tool of choice for studying plant response to temperature variation.

INTRODUCTION
Plant diversity represents a huge reservoir of variations maintained by different evolutionary processes, such as natural or artificial selection. The intraspecific natural variation, i.e., the within-species genetic and phenotypic variations, may reflect species local adaptations to different natural environments (Clauss and Mitchell-Olds, 2006;Fournier-Level et al., 2011;Hancock et al., 2011). Analysing natural variation in wild species could help to understand plant adaptation to specific natural environments from an ecological and evolutionary point of view (Mitchell-Olds and Schmitt, 2006).
For instance, with altitudinal gradients, plants must cope with multiple environmental variations including the decrease of temperatures and air humidity, the increase of UV radiation and a diminution of atmospheric pressure with rising altitudes (Körner, 2007). In response to such climatic variations, plants must tightly regulate their physiological processes and modify their phenotypic traits. Nevertheless, plants need to deal with the trade-off between phenotypic variations necessary to survive and the cost inherent to these changes (Murren et al., 2015). In the context of global climatic changes (altered seasons, temperature changes, occurrence of freezing stress), studying phenotypic variations along altitudinal gradients will help to understand plant ability to cope, acclimate and adapt to these changes (Warren, 1998). Monitoring phenotypic variations following multi-climatic changes is challenging, but plants response to temperature changes can easily be assayed across populations. Seed germination, which is crucial to assure the next generation, is dependent of environmental variations (Penfield and MacGregor, 2017). The rosette and stem growth parameters as well as chlorophyll and anthocyanin contents can be used to illustrate plant responses (Duruflé et al., 2017).
Arabidopsis thaliana is an annual plant species distributed worldwide. It is an interesting model for plant molecular sciences, as well as for ecological and evolutionary studies. Several studies about natural genetic variation of A. thaliana have been published at European distribution scale. A large genetic diversity of A. thaliana with a strong geographic structure was observed and brought the hypothesis of multiple Iberian glacial refuges that contributed to the post-glacial recolonization of Europe by plants (Beck et al., 2008;1001Genomes Consortium, 2016. At a local scale, genetic and phenotypic trait variations were observed across altitudinal gradients between 15 populations from northeast of Spain (Botto, 2015), 10 populations from Norway (Lewandowska-Sabat et al., 2017) and 5 populations from north Italian alps (Günther et al., 2016).
The Pyrenean mountains being a physical barrier between the Iberian Peninsula and the north of Europe have probably been crucial in A. thaliana genetic structure (Lee et al., 2017). In addition, the 30 independent Pyrenean populations collected from separated and deep mountain valleys could reflect local environment adaptive divergence related with genetic variations. As the impact of contrasted temperatures on A. thaliana growth has been addressed in a limited number of studies, the Pyrenean populations, Col and Sha have been grown at 22 • C (optimal growth condition for Col) and 15 • C (potentially optimal for some Pyrenean populations based on WorldClim data). Contrasting temperature regimes using 15 • C as cool condition have been recently tested and have demonstrated the leaves plasticity of two A. thaliana ecotypes Stewart et al., 2016).
The class III peroxidases (CIII Prx) have only been detected in green plants, they constitute a large multigenic family in land plants and due to their cell wall location, they could have important role during cell wall dynamic (Passardi et al., 2004). In addition, CIII Prxs sequences are sufficiently conserved between species; hence, they can be used as genetic markers to establish the structure of populations (Gulsen et al., 2010;Nemli et al., 2014;Pinar et al., 2016). Furthermore, CIII Prxs have already been used to study evolutionary relationships and genotypic diversity on an intra and inter-specific basis (Gulsen et al., 2010;Uzun et al., 2014).
Another valuable genetic marker is ribosomal RNA (rRNA) genes or rDNA. In plants, two major classes of nuclear rDNA can be distinguished: 45S rDNA and 5S rDNA (Layat et al., 2012). The genome of A. thaliana contains hundreds of 45S rDNA units, each of them encoding the 18S, 5.8S and the 25S rRNA genes separated by external (5 ETS and 3 ETS) and internal transcribed spacer (ITS1 and ITS2). The characterization of the 45S rDNA 3 ETS sequence from various A. thaliana populations originating from different geographic regions has revealed the existence of one to four 45S rDNA variants (Pontvianne et al., 2010;Abou-Ellail et al., 2011;Chandrasekhara et al., 2016). Besides, rDNA might have major impact on the whole genome organization in natural populations (Long et al., 2013).
In the present study, (i) the genetic diversity and structure of A. thaliana populations from the French part of Pyrenees mountains have been tested using 1001 genome accessions as reference, (ii) the phenotypic variability in these populations at 22 and 15 • C has been analyzed and (iii) the ability of a given population to acclimate or adapt to different growth temperatures has been evaluated via integrative approaches. This study highlights strong variations in the phenotypic response to environmental constraints, notably temperature variation, which seems to be decoupled from population structure of A. thaliana in the Pyrenees mountains.

Population Sampling, Growth Conditions and Seed Availability
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 N: Number of individuals per population; Localization (GPS coordinates) and environmental variables of the populations of A. thaliana analyzed in this study. Geolocalisation and altitude were determined during the sampling. We used the average of the monthly data to calculate the annual minimum, mean and maximum temperatures and annual precipitation obtained from the WorldClim project (Hijmans et al., 2005; www.worldclim.org). Annual UV radiation was obtained from the Photovoltaic Geographical Information System of the European Community (Huld et al., 2012). Climate PC1 values are the position of populations on the first principal component of the PCA of bioclimatic variables that allows to classify populations in function of their environmental conditions (Supplementary Figure S1). 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 R 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 Urdo 100 3.01 1/1/1/1/1/1/1/1/1/1/1/1 975-999 The following parameters are shown for the populations: Hd, haplotype diversity (Hd, Nei, 1987): the probability that two randomly chosen haplotypes are different in the sample; SNP (%): percentage of SNP found on the concatenated sequences (including the 5 CIII Prx loci) per population; haplotypic frequency distribution: number of individuals belonging to the same haplotype within a population (related to Hd); identity (%): minimum and maximum percentage of identity within population of the concatenated sequences graphically represented in Figure 2.
5 plants per replicate. Plants were grouped by population but their positions were randomized between the three replicates and during the experiments.

Climatic Data
Climatic variables were obtained from WorldClim dataset 1 (Hijmans et al., 2005). The values used are the mean of 30 years (1960 to 1990) with a resolution of ca. 1 km 2 per grid cell. Solar radiation reading were obtained from the Photovoltaic Geographical Information System of the European Communities (Huld et al., 2012)

DNA Extraction
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 3 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.

Genetic Analysis
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 6 iterations and a 10 5 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).

Phylogenetic Analysis by Maximum Likelihood Method
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).

Germination Assays
Seeds were sterilized [30% bleach (v/v), 0.1% triton 100× (v/v)] and sown in Petri dishes containing half-concentrated MS medium (Murashige and Skoog, 1962) supplemented with 1% agar. After 48 h of stratification at 4 • C in darkness, seeds were incubated in a germination cabinet at 22 and 15 • C, 100% of humidity and 40 µmol photons m −2 s −1 continuous white light. The Petri dishes were observed using a Zeiss Axiozoom V16 stereomicroscope. The seed pictures were analyzed using ImageJ. Testa and endosperm ruptures (TR and ER) were quantified, using the testa integrity and the radicle tip protrusion as markers, respectively (Supplementary Figure S3A; Lariguet et al., 2013). The percentages of TR and ER were obtained by counting at least 100 seeds from three independent biological replicates. The results were expressed as the mean percentage of the total seed number for each independent experiment.

Rosettes and Stems Phenotyping
Rosettes grown at 22 and 15 • C were analyzed at 4-and 6-weeks, respectively. These durations correspond to the necessary times for the two Col and Sha accessions to reach the bolting stage at these temperatures (Duruflé et al., 2017). Rosette diameter and fresh weight were determined and the number of leaves was counted. The time to reach bolting stage was checked for each population/accession in a second experiment. From this point, stems length was measured up to the final height. When at least two consecutive stem measurements are equivalent, stem growth was considered as over and stem diameter and length was measured and the lateral stems/cauline leaves were counted. To identify the maximum growth speed of the floral stems, a logistic function was used (Krislov Morris and Kuhn Silk, 1992;Lièvre et al., 2016). The data-fitting curve was created for each stem length kinetic. The logistic function used to model stem growth is given by [Y = Asym/(1 + exp ((xmid -log(t))/scal))] and is parameterized by three parameters corresponding to the final stem length (Asym), the time corresponding to the inflexion point (xmid) and the characteristic growth speed (scal). The estimated parameter was used to characterize individual stem growth speed.

Chlorophyll and Anthocyanin Content
Whole rosettes were ground in liquid nitrogen and used for anthocyanin and chlorophyll extraction. Ground material (0.1 g) was vigorously vortexed with 1 mL of 80% acetone solution for 5 min or with 1 mL of 95% ethanol/1% HCl and stored at 5 • C for 24 h in darkness for chlorophyll or anthocyanin contents, respectively. The samples were centrifuged for 10 min at 1000 g. The absorbance of chlorophyll α and β was measured in the supernatant, respectively at 663 and 647 nm. Chlorophyll concentration was calculated using the following formula: (A 663 × 7.15 + A 647 × 18.71)/mg of fresh material = µg chlorophyll/mg fresh material (Cosio and Dunand, 2010). The anthocyanin content was measured in the supernatant at 530 and 657 nm. Anthocyanin concentration was calculated using the following formula: (A 530 -0.25 × A 657 )/mg of fresh material) (Duruflé et al., 2017).

Statistical Analysis
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 4 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 5 . 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 6 .

Local Environmental Conditions Are Highly Contrasted in the Pyrenees Mountains
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 2 ).
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.

Genetic Structure Analysis Suggests a Specific Lineage in the Pyrenees Mountains
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 outgroups (10 individuals of Lant, and one Col and Sha) and those of 22 accessions available from the 1001 Genomes Project (1001Genomes 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 (nonadmixed) 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). In order to confirm that these 9 populations were homogeneous (Camu, Gava, Grip, Hern, Hosp, Jaco, Pont, Prad and Roch), the size variation of the 3 ETS 45S rRNA gene was used. Four classes or variants exist in Col; three of these variants are abundant (VAR1∼50%; VAR2; VAR3), and one is relatively rare (VAR4) (Pontvianne et al., 2010;Abou-Ellail et al., 2011). However, the majority of the accessions display only one (VAR3 in Sha for example) or two variants (Hinoux and Saez, unpublished observations). Figure 3 shows the different 45S rDNA (R) genotypes found in the 9 analyzed populations. VAR4 was not taken into account to define the R genotypes due to its variability and low frequency. For example, we distinguished R5.1, R5.2 and R5.3 dependant on the relative abundance of VAR1 and VAR3. Pie charts illustrate the proportion of these R genotypes for each population and confirm that the 9 populations (Camu, Gava, Grip, Hern, Hosp, Jaco, Pont, Prad and Roch) are effectively 90-100% homogeneous compared to the out-group Lant population using both ETS and CIII Prxs markers.

The Germination Speed Varies Between Pyrenean Populations
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.
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.

Cold Growth Conditions Induced Various Morphological and Physiological Responses
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 (Figures 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).

Multilevel Analyses Reveal Phenotypic Variations in Response to Growth Temperature
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 (Figures 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.

Effect of Genetic and Environmental Data on Phenotypic Traits
The effect of genetic clusters or environmental data on phenotypic trait variation was evaluated by Kruskal-Wallis non-parametric tests Figures 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).

Contrasted Climatic Conditions for the Pyrenees Populations
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. Argu.

Specific Genetic Lineages Observed in the Pyrenees Mountains
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  Figure S2C) and colored as STRUCTURE analysis (Figure 1) with the admixed populations as genetic cluster 3 in gray. In the column 3, the populations are classified according to their environmental conditions considering as plain, medium and high altitudes (Climate PC1 rank). (B) Results of the Kruskal-Wallis nonparametric test between variables and classification according to the genetic and the environmental groups. Color code represents the p-value of the test (Light blue: p ≤ 0.1; Dark blue: p ≤ 0.5). (C) Stem diameter per genetic cluster and for the two temperatures. (D) Rosette diameter per genetic cluster and for the two temperatures. (E) Chlorophyll content in rosettes per genetic cluster and for the two temperatures. Mean ( ± SEM). p-values (Wilcoxon rank sum test): " * ", p ≤ 0.05; ".", p ≤ 0.1. TR_24 h: Testa rupture at 24 h; ER_30h: Endosperm rupture after 30 h; R_diam: rosette diameter; R_weight: rosette fresh weight; Nb_leaves: rosette leaf number; R_Chloro: rosette chlorophyll content; R_antho: rosette anthocyanin content; Day_bolting: time to reach the bolting stage; St_cauline: cauline leaf number; St_max: stem length; St_diam: floral stem diameter; St_speed: stem elongation speed; Plasticity: Distance separating each Pyrenean population according to the temperature in the multilevel scaled PCA; Cluster_1: Assignment at the genetic cluster 1 determined by STRUCTURE.
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.

Adaptation of the Germination Speed to Different Climatic Constraints
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.

Contrasted Phenotype Induced in Response to Low Temperature
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 multiclimatic 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.

AUTHOR CONTRIBUTIONS
HD, PR, and CD planned, designed the research and performed the bioinformatics analysis. HSC and MaB performed genetic analyses. HD and DLMM performed the phenotyping analyses. AE and SV performed germination assays. NE and MoB collected of natural population and supervised the genetic analyses. PR performed the DNA extraction and the sequencing data. VD-H, JS-V, and J-PR supervised the rRNA genetic analyses and analyzed the data. HD and SD designed the statistical analyses and analyzed the data. HD and CD wrote the article. All the authors read and approved the final manuscript.

FUNDING
This work was supported by the French Laboratory of Excellence project "TULIP" (ANR-10-LABX-41) and IdEx (ANR-11-IDEX-0002-02). HD was supported by an Occitanie Region fellowship and the Federal University of Toulouse.

ACKNOWLEDGMENTS
We are thankful to the Paul Sabatier-Toulouse 3 University and to the Centre National de la Recherche Scientifique (CNRS) for granting their work. We are also grateful to the National Botanical Conservatory of Pyrenees and the Midi-Pyrénées Region (CBNPMP) and to the National Park of the Pyrenees (PNP) for providing us with the population coordinates and giving us sampling authorization within the park limitations. Thanks to Etienne Florence and Alan Griffaut (PNP) for their help during the sampling of populations, to Nicolas Pouilly (LIPM, Toulouse) for advice and assistance with the large scale extraction of Arabidopsis DNA, to Catherine Mathé (LRSV, Toulouse) for her help in choosing the appropriate peroxidase gene markers, to Sarah Pignoly (LRSV, Toulouse) for initial sequences management, Fabrice Roux and Léa Frachon (LIPM, Toulouse) for providing the Lant population and Sophie Brando and Marie Morlans (LGDP, Perpignan) for their technical assistance.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.00430/ full#supplementary-material FIGURE S1 | (A) Primers set used in this study and (B) positions of the amplified sequence on the genes. Gray boxes represent exons, empty boxes introns, and bold lines 5 UTR.
FIGURE S2 | Detail of the steps for the graphical method allowing detection of the number of clusters K estimated with STRUCTURE software. (A) Mean L (K) (±SD) over 10 runs for each K value and (B) delta K calculated as described (Evanno et al., 2005). The modal value of this distribution is the optimal K or the uppermost level of structure, here K = 2. (C) Cluster assignment of Pyrenees populations based on majority membership proportion inferred from K = 2. Each population is represented by a single vertical bar divided into colored segments (green and pink) representing the estimated membership proportions of genetic clusters (K) fitted in the model. Dashed lines corresponded to the delimitation of membership with the 0.8 threshold. Independent batches have been analyzed (mean ± SD) and significant differences between temperature data are shown with Student test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; "·", p ≤ 0.1.  Measurements of the number of leaves per rosettes were done on the day of the harvest (4 and 6 weeks for 22 and 15 • C, respectively). Populations were classified according to the data gathered for plants grown at 22 • C. Three independent batches have been analyzed (mean ± SD) and significant differences between temperature data are shown with Student's test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; ".", p ≤ 0.1). Duncan's test were performed for each growth condition: (B) 22 • C and (C) 15 • C. Different letters indicate significant differences between populations with p ≤ 0.01. Measurements of the weight of rosettes were done on the day of the harvest (4 and 6 weeks for 22 and 15 • C, respectively). Populations were classified according to the data gathered for plants grown at 22 • C. Three independent batches have been analyzed (mean ± SD) and significant differences between temperature data are shown with Student's test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; ".", p ≤ 0.1). Duncan's test were performed for each growth condition: (B) 22 • C and (C) 15 • C. Different letters indicate significant differences between populations with p ≤ 0.01. FIGURE S8 | The diameter of rosettes increases at 15 • C as compared to 22 • C. (A) Measurements of the rosette diameter were done on the day of the harvest (4 and 6 weeks for 22 and 15 • C, respectively). Populations were classified according to the data gathered for plants grown at 22 • C. Three independent batches were analyzed (mean ± SD) and significant differences between temperature data are shown with Student's test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; ".", p ≤ 0.1). Duncan's test were performed for each growth condition: (B) 22 • C and (C) 15 • C. Different letters indicate significant differences between populations with p ≤ 0.01. FIGURE S9 | The time to growth at the bolting stage increases at 15 • C as compared to 22 • C. (A) Different timing to reach to the bolting stage of plant grown in contrasted temperature conditions (mean ± SD) and significant differences between temperature data are shown with Student's test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; ".", p ≤ 0.1). Populations were classified according to the data gathered for plants grown at 22 • C. (B) Comparison between the time to growth at the bolting stage for population grown at 22 and 15 • C. Duncan's test were performed for each growth condition: (C) 22 • C and (D) 15 • C. Different letters indicate significant differences between populations with p ≤ 0.01. Three independent batches have been analyzed.
FIGURE S10 | The anthocyanin content increases in rosettes at 15 • C as compared to 22 • C. (A) Measurements of the anthocyanin content were done on the day of the harvest (4 and 6 weeks for 22 and 15 • C, respectively). Populations were classified according to the data gathered for plants grown at 22 • C. Three independent batches have been analyzed (mean ± SD) and significant differences between temperature data are shown with Student's test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; ".", p ≤ 0.1). Duncan's test were performed for each growth condition: (B) 22 • C and (C) 15 • C. Different letters indicate significant differences between populations with p ≤ 0.01. FIGURE S11 | Chlorophyll content is more constant in rosettes of different population grown at 15 • C as compared to 22 • C. (A) Measurements of the chlorophyll content were done on the day of the harvest (4 and 6 weeks for 22 and 15 • C, respectively). Populations were classified according to the data gathered for plants grown at 22 • C. Three independent batches have been analyzed (mean ± SD) and significant differences between temperature data are shown with Student's test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; ".", p ≤ 0.1). Duncan's test were performed for each growth condition: (B) 22 • C and (C) 15 • C. Different letters indicate significant differences between populations with p ≤ 0.01. (D) Average of chlorophyll content in all the populations at 22 and 15 • C growing conditions. FIGURE S12 | The number of cauline leaves in stems increases at 15 • C as compared to 22 • C. (A) Measurements of the number of cauline leaves were done on the day of the harvest. Populations were classified according to the data gathered for plants grown at 22 • C. Three independent batches have been analyzed (mean ± SD) and significant differences between temperature data are shown with Student's test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; ".", p ≤ 0.1). Duncan's test were performed for each growth condition: (B) 22 • C and (C) 15 • C. Different letters indicate significant differences between populations with p ≤ 0.01. FIGURE S13 | The final stems length increases at 15 • C as compared to 22 • C. (A) Measurements of final stems length were done on the day of the harvest. Populations were classified according to the data gathered for plants grown at 22 • C. Three independent batches have been analyzed (mean ± SD) and significant differences between temperature data are shown with Student's test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; ".", p ≤ 0.1). Duncan's test were performed for each growth condition: (B) 22 • C and (C) 15 • C. Different letters indicate significant differences between populations with p ≤ 0.01. FIGURE S14 | The diameter of stems increases at 15 • C as compared to 22 • C. (A) Measurements of the diameter of stems were done on the day of the harvest. Populations were classified according to the data gathered for plants grown at 22 • C. Three independent batches were analyzed (mean ± SD) and significant differences between temperature data are shown with Student's test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; ".", p ≤ 0.1). Duncan's test were performed for each growth condition: (B) 22 • C and (C) 15 • C. Different letters indicate significant differences between populations with p ≤ 0.01. FIGURE S15 | The stems growth speed increases at 15 • C as compared to 22 • C. (A) Different stems growth speed of plants grown in contrasted temperature conditions and significant differences between temperature data are shown with Student's test: " * * * ", p ≤ 0.001; " * * ", p ≤ 0.01; " * ", p ≤ 0.05; ".", p ≤ 0.1). Populations were classified according to the data gathered for plants grown at 22 • C. (B) Comparison between the stems growth speed at 22 and 15 • C. Duncan's test were performed for each growth condition: (C) 22 • C and (D) 15 • C. Different letters indicate significant differences between populations with p ≤ 0.01. Three independent batches have been analyzed (mean ± SD).