**2. Results**

#### *2.1. Responses of Dormancy Traits of Medicago Accessions to Experimental Temperature Treatments*

Most dormancy traits exhibited a near normal distribution, and a wide range of variability (Figure 1A; Supplementary data Figure S1). Final PY dormancy (FPYD), a proportion of dormant seeds after 88 days of incubation onto water-saturates, ranged from 34% to 100% with mean 80% (SD = 15) at 25/15 ◦C treatment, and from 4% to 94% with mean 60% (SD = 19) at 35/15 ◦C treatment. Comparison of responses of each accession to two temperature treatments showed a remarkable e ffect of larger temperature alternation on dormancy release in a majority of accessions (Figure 1B). The germination pattern (area under curve, AUC M) ranged from 3 to 79, and, similarly to FPYD, larger temperature alternation increased the dormancy release (AUC25: mean ± SD 24 ± 13, range 0–79; AUC35: 34 ± 17, range 5–82), except for some accessions (16%) where the di fferential (AUC35-25) was negative (Figure 1A and Figure S1). Both phenotypic plasticity indexes, based on the minimum and the maximum value among the two temperature treatments divided by the maximum value (PI), showed large ranges with mean PIAUC being slightly higher (0.56 ± 0.29, range 0.00–1.00) than mean PIPY (0.43 ± 0.23, range 0.00–0.91) (Figure 1A and Figure S1). All dormancy traits were moderately to strongly correlated (up to |0.75|, excluding FPYD M and AUC M with some correlations up to |0.94|), except for PIAUC, which was significantly correlated only with AUC35 and AUC25 (Figure 1A; Supplementary data Figure S2, S3).

**Figure 1.** Correlations among dormancy release traits. (**A**) Correlation chart of dormancy release traits and ordination scores of environmental principal component analysis, PCA (first two ordination axes; PCA1, PCA2; see Figure 2A,B). The distribution of each variable is shown on the diagonal. Bellow the diagonal the bivariate scatter plots with a fitted smooth line (loess) are displayed. Above the diagonal the value of the Pearson correlation coefficient plus the significance level as stars are displayed (\* *p* ≤ 0.05, \*\* *p* ≤ 0.01, \*\*\* *p* ≤ 0.001). (**B**) Relationship between final physical dormancy (PY) dormancy of each accession under two temperature treatments (FPYD35 and FPYD25).

**Figure 2.** Principal component analysis (PCA) of selected bioclimatic and soil variables of *Medicago* accessions and multiple correlations of dormancy traits with ordination axes. (**A**,**B**) Principal component analysis (PCA) of selected bioclimatic and soil variables of Medicago accessions. Each accession is classified according to cluster analysis of environmental variables into one of four clusters (see Methods). The ellipses were created based on a model of bivariate normal distribution of the cluster class symbols (estimated from a variance–covariance matrix of their X and Y coordinates) to cover 95% of that distribution' cases. A comparison of selected environmental variables among clusters is shown in Supplementary fileSupplementary file Tables S1 and S2. Vectors of geographic variables (latitude, longitude) were added into the diagram after PCA to visualize spatial gradients of environment. Variables BIO14 and 18 were log(x+1) transformed before analyses. (**C**) Spatial autocorrelation diagram of Moran's I for the first two ordination axes of PCA (PCA1, PCA2). Mean ± 95% CI of I for respective distance class is calculated. (**D**) Multiple correlations of dormancy traits with the first and the second ordination axes of the environmental PCA. Each arrow points in the direction of the steepest increase of the values for corresponding dormancy trait. The angle between arrows indicates the sign of the correlation between the variables. The length of the variable arrows is the multiple correlation of that variable with the ordination axes. Dormancy trait significantly correlated (*p* ≤ 0.05, spatial correlation) with any ordination axis has an asterisk.

#### *2.2. Associations of Environmental Gradients with Dormancy Traits of Medicago Accessions*

Principal component analysis (PCA) of a reduced data set containing 14 climatic and eight soil variables revealed two clear environmental gradients (Figure 2A,B). The first ordination axis explained 30.4% of the total variation and can be interpreted as the gradient of aridity that is tightly correlated with latitude (i.e., the north–south gradient). Climatic variables with the highest positive/negative correlation with the fist ordination axis represent temperatures of the warmest month (BIO5; Pearson's r = 0.61 \*\*\*) and the driest quarter (BIO9, r = 0.71 \*\*\*), isothermality (BIO3, r = 0.59 \*\*\*), precipitation of the driest month (BIO14; r = −0.83 \*\*\*) and precipitation of the warmest quarter (BIO18, r = −0.78 \*\*\*).Concerning soil variables, pH index is positively correlated (r = −0.69 \*\*\*) while soil organic carbon content (ORCDRC, r = −0.75 \*\*\*), available soil water capacity (AWCh1, r = −0.75 \*\*\*), and saturated water content (AWCtS, r = −0.79 \*\*\*) are negatively correlated with the first axis. Latitude (r = −0.77 \*\*\*) but not longitude (r = 0.07) is strongly negatively correlated with the first axis. The second ordination axis explained 17.3% of the total variation and can be interpreted as combined gradient of seasonality and inter-annual variability, with weak geographic (i.e., west–east) trend (latitude: r = −0.01, longitude: r = 0.24 \*\*\*). The most correlated variables with the second axis were precipitation seasonality (BIO15, r = 0.61 \*\*\*) and minimal temperature of the coldest month (BIO6, r = 0.63 \*\*\*), but inter-annual variability of temperature (IV BIO1, r = −0.54 \*\*\*) and precipitation (IV BIO12, r = −0.54 \*\*\*) also had high correlation coefficients. Both synthetic environmental variables (PC1, PC2) were spatially structured as revealed by Moran's I correlogram (both *p* < 0.001), showing positive autocorrelation at short and large distance classes and negative autocorrelation at intermediate distance classes (Figure 2C). Inspection of dormancy trait correlations with ordination axes representing synthetic environmental variables showed that only PIPY was significantly correlated with the first ordination axis (r = 0.16 \*), even after correction for spatial autocorrelation (*p* = 0.032). Other dormancy traits did not show any significant correlation with the first two ordination axes of PCA (Figure 2D; Supplementary data Figure S2). Neither dormancy trait showed any spatial autocorrelation (all Moran's I correlograms had *p* > 0.40, not shown).

Separate analyses of relationships between each dormancy trait and each bioclimatic and soil variable showed that only one dormancy trait (PIPY) was significantly correlated with more environmental variables, while other dormancy traits were either not correlated or showed weak correlations with some environmental variables (Supplementary data Figures S2, S3 and S4, Table S4). Specifically, PIPY was clearly related to the gradient of aridity, i.e., PIPY increases with increasing temperatures and decreasing precipitation and decreasing available soil water capacity (Figure 2; Supplementary data Figure S2). However, there were three climatic variables, i.e., IV BIO1, IV BIO5, and IV BIO10, which showed significant correlations with a majority of dormancy traits (Supplementary data Table S5). Specifically, final PY dormancy (FPYDM, FPYD25, FPYD35) slightly increased with increasing inter-annual variation in temperatures of the warmest quarter (all r = ~0.19 \*).

Four macro-environmental groups of Medicago accessions (Figures 2A and 3, Supplementary data Table S3) differed in slopes of the FPYD across two experimental temperature treatments (Figure 4). Considering each experimental year separately, accessions from arid conditions (clusters K1 and K4, Supplementary data Table S3) consistently showed higher FPYD at 25/15 ◦C and lower at 35/15 ◦C. In contrast, FPYD of accessions from K2 (less arid conditions) did not change significantly in response to different temperature treatments (Figure 4).

#### *2.3. Association Analysis of Dormancy Traits*

In order to identify molecular mechanisms underlying physical dormancy and its adaptability, we performed genome-wide association analyses for all dormancy traits (FPYD25, FPYD35, AUC25, AUC35, AUC35-25, PIPY, PIAUC) and three bioclimatic variables (BIO1, BIO9, BIO12) on 178 accessions. Corresponding Manhattan plots for these analyses are provided in Supplementary data Figure S5. Quantile–quantile (Q-Q) plots confirmed that FarmCPU was a more suitable model to perform association studies (Supplementary data Figure S6). Most significant Quantitative Trait Nucleotides (QTNs) were identified with AUC25, AUC35-25, FPYD25, PIAUC and all three bioclimatic variables. To provide a list of significant QTNs, we defined a threshold of 10−<sup>7</sup> (except for PIPY, where we used a threshold of <sup>10</sup>−4). 136 candidate genes were identified as potential regulators of physical dormancy (Supplementary data Table S6). A large proportion of candidate genes was annotated as involved in synthesis of secondary metabolites, in cell wall modification, and hormone regulation. We performed an over-representation analysis with these 136 candidate genes using a hypergeometric test with Bonferroni correction and this revealed three biological functional Gene Ontology (GO) classes statistically overrepresented (Supplementary data Table S7) and acting as potential regulators of dormancy: response to oxidative stress (GO:0006979), oxidation reduction (GO:0055114), and response

to chemical stimulus (GO:0042221). Candidate genes belonging to these three GO classes are indicated in Table 1.

**Figure 3.** Geographic distribution of studied *Medicago truncatula* accessions classified into four clusters based on climatic and soil conditions, using Ward's minimum-variance linkage of Euclidean distance. Grey dots indicate K1, green K2, light blue K3, and yellow K4 cluster, placed on the background of BIO5 (precipitation in the wet quarter).

**Figure 4.** Reaction norms to changes in experimental temperature (25/15 ◦C, 35/15 ◦C treatments) on final PY dormancy of seeds for K1–K4 macro-environmental clusters in three experimental years (2016, 2017, and 2018). Vertical bars indicate ± SE. Asterisk (\* *p* ≤ 0.05 and \*\* *p* ≤ 0.01) indicates significant differences between temperatures for each cluster.


**Table 1.** List of Quantitative Trait Nucleotides (QTNs) identified by genome-wide association (GWA) analysis of each dormancy trait and belonging to one of three biological function over-represented in our complete list of candidate QTNs. Corresponding chromosome locations and p-values of QTNs are indicated as well as
