# Genome-wide association study identifies favorable SNP alleles and candidate genes for frost tolerance in pea

Sana Beji, Véronique Fontaine, Rosemonde Devaux, Martine Thomas, Sandra Silvia Negro, Nasser Bahrman, Mathieu Siol, Grégoire Aubert, Judith Burstin, Jean-Louis Hilbert, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sana Beji, Véronique Fontaine, Rosemonde Devaux, Martine Thomas, Sandra Silvia Negro, et al. Genome-wide association study identifies favorable SNP alleles and candidate genes for frost tolerance in pea. BMC Genomics, 2020, 21 (1), art.536 (000021p.). 10.1186/s12864-020-06928-w. hal-03270992

HAL Id: hal-03270992

https://hal.inrae.fr/hal-03270992

Submitted on 23 Jan 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

Open Access

# Genome-wide association study identifies favorable SNP alleles and candidate genes for frost tolerance in pea

<figure><img src="image_1.png" /></figure>

Check for
updates

Sana Beji¹*, Véronique Fontaine¹, Rosemonde Devaux², Martine Thomas², Sandra Silvia Negro³, Nasser Bahrman¹, Mathieu Siol⁴, Grégoire Aubert⁴, Judith Burstin⁴, Jean-Louis Hilbert¹, Bruno Delbreil¹ and Isabelle Lejeune-Hénaut¹

## Abstract

**Background:** Frost is a limiting abiotic stress for the winter pea crop (*Pisum sativum* L.) and identifying the genetic determinants of frost tolerance is a major issue to breed varieties for cold northern areas. Quantitative trait loci (QTLs) have previously been detected from bi-parental mapping populations, giving an overview of the genome regions governing this trait. The recent development of high-throughput genotyping tools for pea brings the opportunity to undertake genetic association studies in order to capture a higher allelic diversity within large collections of genetic resources as well as to refine the localization of the causal polymorphisms thanks to the high marker density. In this study, a genome-wide association study (GWAS) was performed using a set of 365 pea accessions. Phenotyping was carried out by scoring frost damages in the field and in controlled conditions. The association mapping collection was also genotyped using an Illumina Infinium® BeadChip, which allowed to collect data for 11,366 single nucleotide polymorphism (SNP) markers.

**Results:** GWAS identified 62 SNPs significantly associated with frost tolerance and distributed over six of the seven pea linkage groups (LGs). These results confirmed 3 QTLs that were already mapped in multiple environments on LG III, V and VI with bi-parental populations. They also allowed to identify one locus, on LG II, which has not been detected yet and two loci, on LGs I and VII, which have formerly been detected in only one environment. Fifty candidate genes corresponding to annotated significant SNPs, or SNPs in strong linkage disequilibrium with the formers, were found to underlie the frost damage (FD)-related loci detected by GWAS. Additionally, the analyses allowed to define favorable haplotypes of markers for the FD-related loci and their corresponding accessions within the association mapping collection.

**Conclusions:** This study led to identify FD-related loci as well as corresponding favorable haplotypes of markers and representative pea accessions that might to be used in winter pea breeding programs. Among the candidate genes highlighted at the identified FD-related loci, the results also encourage further attention to the presence of C-repeat Binding Factors (CBF) as potential genetic determinants of the frost tolerance locus on LG VI.

**Keywords:** Frost damages, Frost tolerance, Genome wide association study (GWAS), Pea (*Pisum sativum* L.), Quantitative trait loci (QTL), Haplotypes of markers, Candidate genes

<footnote>* Correspondence: sana.beji@inrae.fr
¹BioEcoAgro, INRAE, Univ. Liège, Univ. Lille, Univ. Picardie Jules Verne, 2, Chaussée Brunehaut, F-80203 Estrées-Mons, France
Full list of author information is available at the end of the article</footnote>
RESEARCH ARTICLE

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

# Genome-wide association study identifies favorable SNP alleles and candidate genes for frost tolerance in pea

Sana Beji¹*, Véronique Fontaine¹, Rosemonde Devaux², Martine Thomas², Sandra Silvia Negro³, Nasser Bahrman¹, Mathieu Siol⁴, Grégoire Aubert⁴, Judith Burstin⁴, Jean-Louis Hilbert¹, Bruno Delbreil¹ and Isabelle Lejeune-Hénaut¹

## Abstract

**Background:** Frost is a limiting abiotic stress for the winter pea crop (*Pisum sativum* L.) and identifying the genetic determinants of frost tolerance is a major issue to breed varieties for cold northern areas. Quantitative trait loci (QTLs) have previously been detected from bi-parental mapping populations, giving an overview of the genome regions governing this trait. The recent development of high-throughput genotyping tools for pea brings the opportunity to undertake genetic association studies in order to capture a higher allelic diversity within large collections of genetic resources as well as to refine the localization of the causal polymorphisms thanks to the high marker density. In this study, a genome-wide association study (GWAS) was performed using a set of 365 pea accessions. Phenotyping was carried out by scoring frost damages in the field and in controlled conditions. The association mapping collection was also genotyped using an Illumina Infinium® BeadChip, which allowed to collect data for 11,366 single nucleotide polymorphism (SNP) markers.

**Results:** GWAS identified 62 SNPs significantly associated with frost tolerance and distributed over six of the seven pea linkage groups (LGs). These results confirmed 3 QTLs that were already mapped in multiple environments on LG III, V and VI with bi-parental populations. They also allowed to identify one locus, on LG II, which has not been detected yet and two loci, on LGs I and VII, which have formerly been detected in only one environment. Fifty candidate genes corresponding to annotated significant SNPs, or SNPs in strong linkage disequilibrium with the formers, were found to underlie the frost damage (FD)-related loci detected by GWAS. Additionally, the analyses allowed to define favorable haplotypes of markers for the FD-related loci and their corresponding accessions within the association mapping collection.

**Conclusions:** This study led to identify FD-related loci as well as corresponding favorable haplotypes of markers and representative pea accessions that might to be used in winter pea breeding programs. Among the candidate genes highlighted at the identified FD-related loci, the results also encourage further attention to the presence of C-repeat Binding Factors (CBF) as potential genetic determinants of the frost tolerance locus on LG VI.

**Keywords:** Frost damages, Frost tolerance, Genome wide association study (GWAS), Pea (*Pisum sativum* L.), Quantitative trait loci (QTL), Haplotypes of markers, Candidate genes

<footnote>* Correspondence: sana.beji@inrae.fr
¹BioEcoAgro, INRAE, Univ. Liège, Univ. Lille, Univ. Picardie Jules Verne, 2, Chaussée Brunehaut, F-80203 Estrées-Mons, France
Full list of author information is available at the end of the article</footnote>
## Background

In 2018, the world area harvested of pea was ranking behind soybean, common bean, chick pea and cow pea, while the world production of pea was fourth to soybean, common bean and chick pea [1]. Average seed yield worldwide was about 1718 kg/ha in 2018, with the highest yields achieved in the Western European countries [1]. Dry peas are an important nutritional source which provide high quality protein for humans and for animal feeding [2]. In addition to the economic importance of pea seeds, pea crops have beneficial environmental impacts, mainly due to their ability to fix atmospheric nitrogen. They do not need nitrogen fertilizers and therefore help reducing N₂O emissions [3]. For the past decades, spring sowing has been the most common method of cultivation for dry pea. However, the relatively short duration of the development cycle and various stresses such as biotic stresses, mainly *Aphanomyces* root rot and *Ascochyta* blight, as well as abiotic stresses, particularly hydric stress and high temperatures at the end of the development cycle, are at the origin of grain yield losses and variations [4]. Nowadays, winter peas are being developed in order to obtain higher and more stable yields. They are however limited by freezing temperatures during the winter time and the development of winter pea genotypes able to overcome freezing periods is thus desirable.

Mechanisms of tolerance to freezing temperatures have already been reviewed in many plant species. Plants can tolerate freezing temperatures using non-exclusive strategies: freezing escape and cold acclimation. Indeed, plants can escape freezing stress by delaying sensitive phenological stages, particularly floral initiation and flowering, given that frost sensitivity increases after floral initiation [5, 6]. Cold acclimation is the process by which certain plants increase their frost tolerance in response to low non freezing temperatures [7–9]. The *CBF/DREB* (*C*-repeat Binding Factor/Dehydration Responsive Element Binding) transcription factors have an important role in plant cold acclimation. These genes have been isolated first from *Arabidopsis thaliana* and belong to the *AP2/EREBP* (APETALA2/Ethylene-Responsive Element Binding Protein) family of transcription factors [10, 11]. In *Arabidopsis thaliana*, the CBF pathway is characterized by rapid cold induction of *CBF* genes by altering the expression of *CBF*-targeted genes, known as the CBF regulon, which in turn contribute to an increase in freezing tolerance [12]. Many studies have reported the significant role of *CBF* genes in freezing tolerance of herbaceous and woody plant species. Among these studies, the biological role of the CBF pathway for freezing tolerance has also been underlined by the colocalization of *CBF* genes with freezing tolerance quantitative trait loci (QTL) (Arabidopsis: [13], temperate cereals: [14–17], forage grasses: [18], legumes: [19]).

Moreover, within the temperate cereals, *CBF* genes underlying FR2, a major homeologous frost tolerance QTL in barley, diploid and hexaploid wheat, are known to be structured in a cluster of tandemly duplicated genes [20]. In legumes a similar feature was found in *Medicago truncatula*, where a cluster of 12 tandemly duplicated *CBF* genes was shown to match with a major freezing tolerance QTL on chromosome 6 [19].

The identification of genomic regions controlling frost tolerance has initially been completed for cultivated species through the assessment of mapping populations and QTL mapping. In pea, QTL mapping studies for frost damages have been conducted in multiple field environments as well as in controlled conditions [21–23]. QTLs were detected using two populations of recombinant inbred lines (RILs), namely Pop2 and Pop9, respectively derived from crosses between Champagne (frost tolerant) and Térèse (frost sensitive) [21, 22] and China (frost tolerant) and Caméor (frost sensitive) [23]. Four common QTLs were detected within both populations. The corresponding regions were located on linkage groups (LG) III (two independent positions), V and VI. They explained altogether a major part of the phenotypic variation and were repeatable across environmental conditions. Three other loci were specific to one or other of the two populations and detected in fewer environments (Pop2: one locus on LG I, Pop9: one locus on LG VII, detected in one environment; Pop9: one locus on LG V, detected in two environments).

The resolution and accuracy with which QTL mapping can identify causal genetic determinism of the considered traits is limited by the high confidence intervals and the relatively low total number of recombination events in bi-parental mapping populations. In addition to QTL mapping, association studies have emerged as a complementary approach to dissect quantitative traits by exploiting natural genetic diversity and ancestral recombination events present in germplasm collections [24]. Used for more than two decades in human genetic research, association genetics have been adapted for genetic dissection in plants, taking advantage of the development of high-throughput genotyping resources in numerous species [24]. Genome-wide association studies (GWAS) aim at identifying genetic markers strongly associated with quantitative traits by using the linkage disequilibrium (LD) between candidate genes and markers. They rely on high-density genetic maps allowing an increased resolution of detection generating more precise QTL positions than bi-parental QTL mapping and give access to multiple allelic variation through the exploration of natural genetic diversity. In addition, GWAS can offer a powerful genomic tool for breeding plants by the identification of associated markers tightly linked to targeted genomic regions which can be used for marker-assisted selection. In the recent years, GWAS has
## Background

In 2018, the world area harvested of pea was ranking behind soybean, common bean, chick pea and cow pea, while the world production of pea was fourth to soybean, common bean and chick pea [1]. Average seed yield worldwide was about 1718 kg/ha in 2018, with the highest yields achieved in the Western European countries [1]. Dry peas are an important nutritional source which provide high quality protein for humans and for animal feeding [2]. In addition to the economic importance of pea seeds, pea crops have beneficial environmental impacts, mainly due to their ability to fix atmospheric nitrogen. They do not need nitrogen fertilizers and therefore help reducing N₂O emissions [3]. For the past decades, spring sowing has been the most common method of cultivation for dry pea. However, the relatively short duration of the development cycle and various stresses such as biotic stresses, mainly *Aphanomyces* root rot and *Ascochyta* blight, as well as abiotic stresses, particularly hydric stress and high temperatures at the end of the development cycle, are at the origin of grain yield losses and variations [4]. Nowadays, winter peas are being developed in order to obtain higher and more stable yields. They are however limited by freezing temperatures during the winter time and the development of winter pea genotypes able to overcome freezing periods is thus desirable.

Mechanisms of tolerance to freezing temperatures have already been reviewed in many plant species. Plants can tolerate freezing temperatures using non-exclusive strategies: freezing escape and cold acclimation. Indeed, plants can escape freezing stress by delaying sensitive phenological stages, particularly floral initiation and flowering, given that frost sensitivity increases after floral initiation [5, 6]. Cold acclimation is the process by which certain plants increase their frost tolerance in response to low non freezing temperatures [7–9]. The *CBF/DREB* (*C*-repeat Binding Factor/Dehydration Responsive Element Binding) transcription factors have an important role in plant cold acclimation. These genes have been isolated first from *Arabidopsis thaliana* and belong to the *AP2/EREBP* (APETALA2/Ethylene-Responsive Element Binding Protein) family of transcription factors [10, 11]. In *Arabidopsis thaliana*, the CBF pathway is characterized by rapid cold induction of *CBF* genes by altering the expression of *CBF*-targeted genes, known as the CBF regulon, which in turn contribute to an increase in freezing tolerance [12]. Many studies have reported the significant role of *CBF* genes in freezing tolerance of herbaceous and woody plant species. Among these studies, the biological role of the CBF pathway for freezing tolerance has also been underlined by the colocalization of *CBF* genes with freezing tolerance quantitative trait loci (QTL) (Arabidopsis: [13], temperate cereals: [14–17], forage grasses: [18], legumes: [19]). Moreover, within the temperate cereals, *CBF* genes underlying FR2, a major homeologous frost tolerance QTL in barley, diploid and hexaploid wheat, are known to be structured in a cluster of tandemly duplicated genes [20]. In legumes a similar feature was found in *Medicago truncatula*, where a cluster of 12 tandemly duplicated *CBF* genes was shown to match with a major freezing tolerance QTL on chromosome 6 [19].

The identification of genomic regions controlling frost tolerance has initially been completed for cultivated species through the assessment of mapping populations and QTL mapping. In pea, QTL mapping studies for frost damages have been conducted in multiple field environments as well as in controlled conditions [21–23]. QTLs were detected using two populations of recombinant inbred lines (RILs), namely Pop2 and Pop9, respectively derived from crosses between Champagne (frost tolerant) and Térèse (frost sensitive) [21, 22] and China (frost tolerant) and Caméor (frost sensitive) [23]. Four common QTLs were detected within both populations. The corresponding regions were located on linkage groups (LG) III (two independent positions), V and VI. They explained altogether a major part of the phenotypic variation and were repeatable across environmental conditions. Three other loci were specific to one or other of the two populations and detected in fewer environments (Pop2: one locus on LG I, Pop9: one locus on LG VII, detected in one environment; Pop9: one locus on LG V, detected in two environments).

The resolution and accuracy with which QTL mapping can identify causal genetic determinism of the considered traits is limited by the high confidence intervals and the relatively low total number of recombination events in bi-parental mapping populations. In addition to QTL mapping, association studies have emerged as a complementary approach to dissect quantitative traits by exploiting natural genetic diversity and ancestral recombination events present in germplasm collections [24]. Used for more than two decades in human genetic research, association genetics have been adapted for genetic dissection in plants, taking advantage of the development of high-throughput genotyping resources in numerous species [24]. Genome-wide association studies (GWAS) aim at identifying genetic markers strongly associated with quantitative traits by using the linkage disequilibrium (LD) between candidate genes and markers. They rely on high-density genetic maps allowing an increased resolution of detection generating more precise QTL positions than bi-parental QTL mapping and give access to multiple allelic variation through the exploration of natural genetic diversity. In addition, GWAS can offer a powerful genomic tool for breeding plants by the identification of associated markers tightly linked to targeted genomic regions which can be used for marker-assisted selection. In the recent years, GWAS has
been conducted in many plant species to dissect complex quantitative traits including winter survival and frost tolerance [25–33]. High throughput genotyping resources now available in pea [34] have also allowed to carry out GWAS in order to dissect the genetic determinism of resistance to *Aphanomyces euteiches*, plant architecture and frost tolerance. Desgroux et al. [35] studied associations for resistance to *A. euteiches* in 175 pea lines using a high-density SNP genotyping consisting of 13.2 K SNPs from the developed GenoPea Infinium® BeadChip [36]. Several markers significantly associated with resistance to *A. euteiches* harboring relevant putative candidate genes were identified. Significantly associated markers also allowed to refine the confidence interval of QTLs previously detected in biparental populations. Using the same SNP resource and a collection of 266 pea accessions, including the 175 former lines, Desgroux et al. also identified genomic intervals significantly associated with plant architecture and resistance to *A. euteiches*, of which 8 were overlapping for both traits [37]. In a different genetic background composed by a set of 672 pea accessions genotyped with 267 simple sequence repeat (SSR) markers, Liu et al. [38] detected 7 SSRs significantly associated with frost tolerance of which one was located on LG VI and was shown to colocalize with a gene involved in the metabolism of glycoproteins in response to chilling stress.

The present study aimed at reconsidering the regions of the genome that control frost tolerance in pea, taking advantage of genome-wide distributed SNPs generated from the 13.2 K GenoPea Infinium® BeadChip [36].

## Results

### Statistical analyses of phenotypic data

In order to undertake genome-wide association analyzes, a collection of 365 pea accessions, here after referred to as the association mapping collection or the collection, was phenotyped for frost damages (FD) under field and controlled conditions. Statistical analyses of frost damages scores showed highly significant genetic variation for all studied traits and the coefficient of genetic variation ranged between 37.6 for FD\_Field\_date4 and 74.4

for FD\_Field\_date2 (Table 1). The estimates of broad sense heritabilities (H²) were high for all traits, varying from 0.84 to 0.89. In addition, frost damages observed in controlled conditions and in the field for the date 3 and 4 showed the highest mean scores of 2.8, 2.4 and 2.4 respectively (Table 1). Frequency distributions of best linear unbiased predictions (BLUPs) values for each trait tended to fit normal curves within the collection (Additional file 1).

### SNP genotyping of the association mapping collection

After quality control, the genotyping data comprised a total of 10,739 polymorphic SNPs with imputed missing data and a minor allele frequency (MAF) greater than 5%. Each linkage group contained 1533 SNPs on average. The distribution of SNPs varied within and between LGs (Additional files 2 and 3). Linkage group III showed the highest number of SNPs (1888 SNPs), while LG I was the least dense (1220 SNPs). Mean MAF in the association mapping collection varied from 0.29 on LGI to 0.31 on LG V and LG VI (Additional file 3). Only 751 of SNPs had a MAF less than or equal to 0.1 (Additional file 4). The distribution of SNP markers across the different LGs was dense and no gap between adjacent SNPs exceeded 1.7 cM, except on LG I and LG V which presented gaps of 2.3 cM (for the interval position 48.4–50.7 cM) and 3.7 cM (for the interval position 0.1–3.7 cM), respectively (Additional file 2). In addition, the map of the 10,739 SNPs used for GWAS showed an average number of 28 SNPs mapped at the same genetic position.

### LD analysis

The distribution of the estimate of the linkage disequilibrium (LD, r²) along to the genetic position for each linkage group as well as for the whole genome is presented in Additional file 5. The r² value rapidly decreased as the genetic distance increased. The LD decay, estimated as the distance for which r² decreases to half of its maximum level (0.22), was equal to 0.9 cM for the whole genome. Considering the LGs individually, the LD decay ranged from 0.3 cM for LG IV to 1.4 cM for LG V.

<table><caption>**Table 1** Statistical parameters of the pea collection for the five observed traits</caption><thead><tr><th>Trait</th><th>Number of accessions</th><th>Number of observations</th><th>Min</th><th>Max</th><th>Mean</th><th>SE</th><th>Vg</th><th>CVg</th><th>H<sup>2</sup></th></tr></thead><tbody><tr><td>FD_CC</td><td>363</td><td>1087</td><td>0.00</td><td>5.00</td><td>2.80</td><td>0.03</td><td>1.25</td><td>39.93</td><td>0.89</td></tr><tr><td>FD_Field_date1</td><td>363</td><td>1089</td><td>0.00</td><td>5.00</td><td>1.43</td><td>0.03</td><td>0.87</td><td>65.23</td><td>0.89</td></tr><tr><td>FD_Field_date2</td><td>363</td><td>1048</td><td>0.00</td><td>5.00</td><td>1.39</td><td>0.03</td><td>1.07</td><td>74.42</td><td>0.89</td></tr><tr><td>FD_Field_date3</td><td>363</td><td>1083</td><td>0.00</td><td>5.00</td><td>2.40</td><td>0.04</td><td>1.65</td><td>53.52</td><td>0.87</td></tr><tr><td>FD_Field_date4</td><td>363</td><td>849</td><td>0.00</td><td>5.00</td><td>2.41</td><td>0.03</td><td>0.82</td><td>37.57</td><td>0.84</td></tr></tbody></table>

Traits are abbreviated as follows: FD\_CC: Frost damages in the controlled conditions experiment, FD\_Field\_date1, FD\_Field\_date2, FD\_Field\_date3 and FD\_Field\_date4: Frost damages in the field experiment at the first, second, third and fourth date of observation, respectively. Minimum (Min), maximum (Max) and mean values, standard error (SE), genetic variance (Vg), coefficient of variation of the genetic variance (CVg) expressed as √(Vg)/mean and broad-sense heritability (H²) are shown for each trait
been conducted in many plant species to dissect complex quantitative traits including winter survival and frost tolerance [25–33]. High throughput genotyping resources now available in pea [34] have also allowed to carry out GWAS in order to dissect the genetic determinism of resistance to *Aphanomyces euteiches*, plant architecture and frost tolerance. Desgroux et al. [35] studied associations for resistance to *A. euteiches* in 175 pea lines using a high-density SNP genotyping consisting of 13.2 K SNPs from the developed GenoPea Infinium® BeadChip [36]. Several markers significantly associated with resistance to *A. euteiches* harboring relevant putative candidate genes were identified. Significantly associated markers also allowed to refine the confidence interval of QTLs previously detected in biparental populations. Using the same SNP resource and a collection of 266 pea accessions, including the 175 former lines, Desgroux et al. also identified genomic intervals significantly associated with plant architecture and resistance to *A. euteiches*, of which 8 were overlapping for both traits [37]. In a different genetic background composed by a set of 672 pea accessions genotyped with 267 simple sequence repeat (SSR) markers, Liu et al. [38] detected 7 SSRs significantly associated with frost tolerance of which one was located on LG VI and was shown to colocalize with a gene involved in the metabolism of glycoproteins in response to chilling stress.

The present study aimed at reconsidering the regions of the genome that control frost tolerance in pea, taking advantage of genome-wide distributed SNPs generated from the 13.2 K GenoPea Infinium® BeadChip [36].

# Results

## Statistical analyses of phenotypic data

In order to undertake genome-wide association analyzes, a collection of 365 pea accessions, here after referred to as the association mapping collection or the collection, was phenotyped for frost damages (FD) under field and controlled conditions. Statistical analyses of frost damages scores showed highly significant genetic variation for all studied traits and the coefficient of genetic variation ranged between 37.6 for FD\_Field\_date4 and 74.4

for FD\_Field\_date2 (Table 1). The estimates of broad sense heritabilities (H²) were high for all traits, varying from 0.84 to 0.89. In addition, frost damages observed in controlled conditions and in the field for the date 3 and 4 showed the highest mean scores of 2.8, 2.4 and 2.4 respectively (Table 1). Frequency distributions of best linear unbiased predictions (BLUPs) values for each trait tended to fit normal curves within the collection (Additional file 1).

## SNP genotyping of the association mapping collection

After quality control, the genotyping data comprised a total of 10,739 polymorphic SNPs with imputed missing data and a minor allele frequency (MAF) greater than 5%. Each linkage group contained 1533 SNPs on average. The distribution of SNPs varied within and between LGs (Additional files 2 and 3). Linkage group III showed the highest number of SNPs (1888 SNPs), while LG I was the least dense (1220 SNPs). Mean MAF in the association mapping collection varied from 0.29 on LGI to 0.31 on LG V and LG VI (Additional file 3). Only 751 of SNPs had a MAF less than or equal to 0.1 (Additional file 4). The distribution of SNP markers across the different LGs was dense and no gap between adjacent SNPs exceeded 1.7 cM, except on LG I and LG V which presented gaps of 2.3 cM (for the interval position 48.4–50.7 cM) and 3.7 cM (for the interval position 0.1–3.7 cM), respectively (Additional file 2). In addition, the map of the 10,739 SNPs used for GWAS showed an average number of 28 SNPs mapped at the same genetic position.

## LD analysis

The distribution of the estimate of the linkage disequilibrium (LD, r²) along to the genetic position for each linkage group as well as for the whole genome is presented in Additional file 5. The r² value rapidly decreased as the genetic distance increased. The LD decay, estimated as the distance for which r² decreases to half of its maximum level (0.22), was equal to 0.9 cM for the whole genome. Considering the LGs individually, the LD decay ranged from 0.3 cM for LG IV to 1.4 cM for LG V.

<table><caption>**Table 1** Statistical parameters of the pea collection for the five observed traits</caption><thead><tr><th>Trait</th><th>Number of accessions</th><th>Number of observations</th><th>Min</th><th>Max</th><th>Mean</th><th>SE</th><th>Vg</th><th>CVg</th><th>H<sup>2</sup></th></tr></thead><tbody><tr><td>FD_CC</td><td>363</td><td>1087</td><td>0.00</td><td>5.00</td><td>2.80</td><td>0.03</td><td>1.25</td><td>39.93</td><td>0.89</td></tr><tr><td>FD_Field_date1</td><td>363</td><td>1089</td><td>0.00</td><td>5.00</td><td>1.43</td><td>0.03</td><td>0.87</td><td>65.23</td><td>0.89</td></tr><tr><td>FD_Field_date2</td><td>363</td><td>1048</td><td>0.00</td><td>5.00</td><td>1.39</td><td>0.03</td><td>1.07</td><td>74.42</td><td>0.89</td></tr><tr><td>FD_Field_date3</td><td>363</td><td>1083</td><td>0.00</td><td>5.00</td><td>2.40</td><td>0.04</td><td>1.65</td><td>53.52</td><td>0.87</td></tr><tr><td>FD_Field_date4</td><td>363</td><td>849</td><td>0.00</td><td>5.00</td><td>2.41</td><td>0.03</td><td>0.82</td><td>37.57</td><td>0.84</td></tr></tbody></table>

Traits are abbreviated as follows: FD\_CC: Frost damages in the controlled conditions experiment, FD\_Field\_date1, FD\_Field\_date2, FD\_Field\_date3 and FD\_Field\_date4: Frost damages in the field experiment at the first, second, third and fourth date of observation, respectively. Minimum (Min), maximum (Max) and mean values, standard error (SE), genetic variance (Vg), coefficient of variation of the genetic variance (CVg) expressed as $\sqrt{Vg}/\text{mean}$ and broad-sense heritability (H²) are shown for each trait
Population structure and kinship analyses

To avoid false positive results in association analysis,
structure and kinship of the association mapping
collection were analyzed using 2962 non-redundant
positions markers. The collection structure was studied
using the discriminant analysis of principal components
(DAPC) method. Following the analysis of the Bayesian

Information Criterion (BIC) profile and using the ‘a-
score’ criterion, the optimal number of clusters was
fixed to 7 (Fig. 1a) and the optimal number of principal
components (PCs) was set to 6 (Fig. 1b). Therefore,
these 6 PCs and 7 clusters were used for discriminant
analysis of principal components. The distribution of
individuals into the 7 clusters is represented along the

<figure><img src="image_2.png" /></figure>
Population structure and kinship analyses

To avoid false positive results in association analysis,
structure and kinship of the association mapping
collection were analyzed using 2962 non-redundant
positions markers. The collection structure was studied
using the discriminant analysis of principal components
(DAPC) method. Following the analysis of the Bayesian

Information Criterion (BIC) profile and using the ‘a-
score’ criterion, the optimal number of clusters was
fixed to 7 (Fig. 1a) and the optimal number of principal
components (PCs) was set to 6 (Fig. 1b). Therefore,
these 6 PCs and 7 clusters were used for discriminant
analysis of principal components. The distribution of
individuals into the 7 clusters is represented along the

<figure><img src="image_4.png" /></figure>
two first axes of DAPC (Fig. 1c). The main passport data (Additional file 6) seemed to be related to the discrimination of the clusters. The cluster 1, comprising mainly wild peas and landraces from Africa and Middle or Far-East, was totally separated from the other clusters (Fig. 1d). The majority of accessions from clusters 2, 5, 6 and 4 were registered as spring sowing types while the winter sowing types were essentially gathered in clusters 3 and 7. These last two clusters differed for their end-use, the cluster 3 being mainly composed of field peas (81%) and the cluster 7 of fodder peas (80%).

The dendrogram using Nei genetic distance among accessions of the association mapping collection revealed also the presence of seven clusters (Additional file 7), as the DAPC method. For 79.61% of the accessions, the assignment to clusters performed by the dendrogram corresponded to the allocation made by the DAPC analysis (Additional file 6).

Within the kinship matrix (K), estimated for the whole genome, 85.7% of the kinship coefficient values were less than 0.1. 1.6% of these values were larger than 0.5. For the seven kinship matrices specific to each linkage group (K<sub>LG</sub>), the kinship coefficient values ranged similarly than for the K matrix (Additional file 8). These results indicated a weak relatedness between accessions and suggested that the majority of the accessions are genetically diverse, which was beneficial for subsequent GWAS mapping. From these results, the two first coordinates of DAPC results (Q matrix) and relatedness matrices (K matrix and K<sub>LG</sub> matrix) were used as covariates for subsequent association analyses.

## Genome-wide association mapping

The comparison of BIC values of the four GWA-models tested, showed that the linear mixed model which included both Q and K<sub>LG</sub> matrices as covariates was the optimal model for the following traits: FD_CC, FD_field_date1, FD_field_date2 and FD_field_date3. Whereas, the best fitting model for FD_field_date4 was the linear mixed model which included only K<sub>LG</sub> matrices (Table 2). The Manhattan and

their corresponding quantile-quantile plots of the association mapping results, run with the best model for each trait, are presented in Figs. 2 and 3. A total of 62 markers were significantly associated with any of the studied traits at the Bonferroni threshold -log10 (p) > 5.33. Frost Damage (FD)-associated markers were distributed on all linkage groups except LGIV. These SNPs exhibited a minor allele frequency ranging from 0.13 to 0.5. The number of markers associated with FD_CC, FD_field_date1, FD_field_date2, FD_field_date3 and FD_field_date4 were 40, 4, 6, 3 and 17, respectively. The highest p values were showed by the significant loci located respectively on LGV (9.67E-11) and LGVI (1.02E-10) (Table 3).

## Favorable alleles exploration for frost tolerance in pea

Twelve LD blocks were defined around the 62 significant FD-associated markers which included all markers in LD ($r^2$ > 0.8) with the FD-associated markers (Additional file 9). FD-associated markers which were not in significant LD with any other SNP, and thus did not constitute a LD block, were also kept for further analysis. Finally, 75 SNPs, covering six FD-related loci distributed on LG I, II, III, V, VI and VII, were kept to identify favorable and unfavorable haplotypes for frost tolerance. For each of the six FD-related loci, marker haplotypes (two to nine) and corresponding representative accessions were identified (Additional file 10). For each locus, the effect of the different allele combinations was tested thanks to a variance analysis and a multiple comparison test of phenotypic mean effects (Additional file 10). These analyses identified 7 favorable haplotypes over the 6 FD-related regions carrying favorable alleles at each FD-associated marker except the haplotypes V.3 and VII.4 which contained each unfavorable allele, among 1 and 3 FD-associated markers respectively. Accessions carrying favorable haplotypes presented lower values of frost damages ranging between 0.00 ± 0.00 (haplotype VII.4 for the trait FD_field_date2) and 2.65 ± 0.05 (haplotype III.1 for the trait FD_CC). Six groups of accessions carrying unfavorable haplotypes were also identified for which 100% of unfavorable alleles were observed over the significant FD-associated markers detected by GWAS.

<table><caption>**Table 2** BIC-based comparison of the four models used to control the rate of false positive associations</caption><thead><tr><th rowspan="2">Trait</th><th colspan="2">LMM1</th><th colspan="2">LMM2</th><th colspan="2">LMM3</th><th colspan="2">LMM4</th></tr><tr><th>-2ln(L)</th><th>BIC</th><th>-2ln(L)</th><th>BIC</th><th>-2ln(L)</th><th>BIC</th><th>-2ln(L)</th><th>BIC</th></tr></thead><tbody><tr><td>FD_CC</td><td>704.71</td><td>4073.94</td><td>701.16</td><td>4070.4</td><td>695.27</td><td>4064.51</td><td>691.68</td><td><strong>4060.92</strong></td></tr><tr><td>FD_field_date1</td><td>658.87</td><td>4028.11</td><td>658.4</td><td>4027.63</td><td>675.69</td><td>4026.92</td><td>656.14</td><td><strong>4025.38</strong></td></tr><tr><td>FD_field_date2</td><td>735.42</td><td>4076.81</td><td>732.06</td><td>4073.45</td><td>731.26</td><td>4072.65</td><td>728.85</td><td><strong>4070.24</strong></td></tr><tr><td>FD_field_date3</td><td>837.83</td><td>4207.06</td><td>833.87</td><td>4203.1</td><td>836.8</td><td>4206.04</td><td>831.51</td><td><strong>4200.74</strong></td></tr><tr><td>FD_field_date4</td><td>482.59</td><td>3489.84</td><td>481.45</td><td><strong>3488.7</strong></td><td>483.23</td><td>3490.48</td><td>482.53</td><td>3489.78</td></tr></tbody></table>

LMM1 linear mixed model including the K kinship matrix, LMM2 linear mixed model including the K<sub>LG</sub> kinship matrix, LMM3 linear mixed model including the K kinship matrix and the population structure matrix (Q), LMM4 linear mixed model including the K<sub>LG</sub> kinship matrix and the population structure matrix (Q). The model with the lowest BIC value (bold) is considered to be the best choice for the target trait
two first axes of DAPC (Fig. 1c). The main passport data (Additional file 6) seemed to be related to the discrimination of the clusters. The cluster 1, comprising mainly wild peas and landraces from Africa and Middle or Far-East, was totally separated from the other clusters (Fig. 1d). The majority of accessions from clusters 2, 5, 6 and 4 were registered as spring sowing types while the winter sowing types were essentially gathered in clusters 3 and 7. These last two clusters differed for their end-use, the cluster 3 being mainly composed of field peas (81%) and the cluster 7 of fodder peas (80%).

The dendrogram using Nei genetic distance among accessions of the association mapping collection revealed also the presence of seven clusters (Additional file 7), as the DAPC method. For 79.61% of the accessions, the assignment to clusters performed by the dendrogram corresponded to the allocation made by the DAPC analysis (Additional file 6).

Within the kinship matrix (K), estimated for the whole genome, 85.7% of the kinship coefficient values were less than 0.1. 1.6% of these values were larger than 0.5. For the seven kinship matrices specific to each linkage group (K<sub>LG</sub>), the kinship coefficient values ranged similarly than for the K matrix (Additional file 8). These results indicated a weak relatedness between accessions and suggested that the majority of the accessions are genetically diverse, which was beneficial for subsequent GWAS mapping. From these results, the two first coordinates of DAPC results (Q matrix) and relatedness matrices (K matrix and K<sub>LG</sub> matrix) were used as covariates for subsequent association analyses.

## Genome-wide association mapping

The comparison of BIC values of the four GWA-models tested, showed that the linear mixed model which included both Q and K<sub>LG</sub> matrices as covariates was the optimal model for the following traits: FD_CC, FD_field_date1, FD_field_date2 and FD_field_date3. Whereas, the best fitting model for FD_field_date4 was the linear mixed model which included only K<sub>LG</sub> matrices (Table 2). The Manhattan and

their corresponding quantile-quantile plots of the association mapping results, run with the best model for each trait, are presented in Figs. 2 and 3. A total of 62 markers were significantly associated with any of the studied traits at the Bonferroni threshold -log10 (p) > 5.33. Frost Damage (FD)-associated markers were distributed on all linkage groups except LGIV. These SNPs exhibited a minor allele frequency ranging from 0.13 to 0.5. The number of markers associated with FD_CC, FD_field_date1, FD_field_date2, FD_field_date3 and FD_field_date4 were 40, 4, 6, 3 and 17, respectively. The highest p values were showed by the significant loci located respectively on LGV (9.67E-11) and LGVI (1.02E-10) (Table 3).

## Favorable alleles exploration for frost tolerance in pea

Twelve LD blocks were defined around the 62 significant FD-associated markers which included all markers in LD ($r^2$ > 0.8) with the FD-associated markers (Additional file 9). FD-associated markers which were not in significant LD with any other SNP, and thus did not constitute a LD block, were also kept for further analysis. Finally, 75 SNPs, covering six FD-related loci distributed on LG I, II, III, V, VI and VII, were kept to identify favorable and unfavorable haplotypes for frost tolerance. For each of the six FD-related loci, marker haplotypes (two to nine) and corresponding representative accessions were identified (Additional file 10). For each locus, the effect of the different allele combinations was tested thanks to a variance analysis and a multiple comparison test of phenotypic mean effects (Additional file 10). These analyses identified 7 favorable haplotypes over the 6 FD-related regions carrying favorable alleles at each FD-associated marker except the haplotypes V.3 and VII.4 which contained each unfavorable allele, among 1 and 3 FD-associated markers respectively. Accessions carrying favorable haplotypes presented lower values of frost damages ranging between 0.00 ± 0.00 (haplotype VII.4 for the trait FD_field_date2) and 2.65 ± 0.05 (haplotype III.1 for the trait FD_CC). Six groups of accessions carrying unfavorable haplotypes were also identified for which 100% of unfavorable alleles were observed over the significant FD-associated markers detected by GWAS.

<table><caption>**Table 2** BIC-based comparison of the four models used to control the rate of false positive associations</caption><thead><tr><th rowspan="2">Trait</th><th colspan="2">LMM1</th><th colspan="2">LMM2</th><th colspan="2">LMM3</th><th colspan="2">LMM4</th></tr><tr><th>-2ln(L)</th><th>BIC</th><th>-2ln(L)</th><th>BIC</th><th>-2ln(L)</th><th>BIC</th><th>-2ln(L)</th><th>BIC</th></tr></thead><tbody><tr><td>FD_CC</td><td>704.71</td><td>4073.94</td><td>701.16</td><td>4070.4</td><td>695.27</td><td>4064.51</td><td>691.68</td><td><strong>4060.92</strong></td></tr><tr><td>FD_field_date1</td><td>658.87</td><td>4028.11</td><td>658.4</td><td>4027.63</td><td>675.69</td><td>4026.92</td><td>656.14</td><td><strong>4025.38</strong></td></tr><tr><td>FD_field_date2</td><td>735.42</td><td>4076.81</td><td>732.06</td><td>4073.45</td><td>731.26</td><td>4072.65</td><td>728.85</td><td><strong>4070.24</strong></td></tr><tr><td>FD_field_date3</td><td>837.83</td><td>4207.06</td><td>833.87</td><td>4203.1</td><td>836.8</td><td>4206.04</td><td>831.51</td><td><strong>4200.74</strong></td></tr><tr><td>FD_field_date4</td><td>482.59</td><td>3489.84</td><td>481.45</td><td><strong>3488.7</strong></td><td>483.23</td><td>3490.48</td><td>482.53</td><td>3489.78</td></tr></tbody></table>

LMM1 linear mixed model including the K kinship matrix, LMM2 linear mixed model including the K<sub>LG</sub> kinship matrix, LMM3 linear mixed model including the K kinship matrix and the population structure matrix (Q), LMM4 linear mixed model including the K<sub>LG</sub> kinship matrix and the population structure matrix (Q). The model with the lowest BIC value (bold) is considered to be the best choice for the target trait
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Manhattan plots of markers associated with the five frost damage traits. The plots show the *p*-values (*p*) for association between a phenotypical trait and each tested marker (expressed as the negative decimal logarithm of *p*, y-axis) plotted against the linkage group position of the marker (x-axis). Plots above the red horizontal line indicate the genome-wide significance with the Bonferroni threshold (–log10 (*p*) > 5.33). **a** is the plot for the evaluation of frost damages in the controlled conditions experiment. **b**, **c**, **d** and **e** are the plots for the four evaluations of FD in the field experiment, corresponding to the 4 dates of damages observation</figcaption></figure>

Typical accessions carrying favorable haplotypes and show-
ing a mean score of frost damage ≤1 were mainly winter
fodder peas (e.g. Black seeded, Champagne, Melrose, Blixt
7) (Additional files 6 and 10). While those carrying unfavor-
able haplotypes with a mean score of frost damage ≥4 were
mainly spring garden peas (e.g. Automobile, Caroubel,
Cennia, Ersling) (Additional files 6 and 10). The sequences
of the 75 SNP markers related to frost tolerance trait are
provided in Tayeh et al. [36] (Table S2).

Candidate genes

The projection of the 75 FD-related markers on the pea
genome assembly [39] allowed to define intervals of 2
Mb on all the pea chromosomes, except the chromo-
some 4 (LG IV). Four FD-related markers were assigned
to unanchored scaffolds which were all less than 2 Mb
long: in that case, annotated genes were listed for the
whole scaffold. A particular case was encountered for
the associated marker PsCam036704_21832_970 which
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Manhattan plots of markers associated with the five frost damage traits. The plots show the *p*-values (*p*) for association between a phenotypical trait and each tested marker (expressed as the negative decimal logarithm of *p*, y-axis) plotted against the linkage group position of the marker (x-axis). Plots above the red horizontal line indicate the genome-wide significance with the Bonferroni threshold (–log10 (*p*) > 5.33). **a** is the plot for the evaluation of frost damages in the controlled conditions experiment. **b**, **c**, **d** and **e** are the plots for the four evaluations of FD in the field experiment, corresponding to the 4 dates of damages observation</figcaption></figure>

Typical accessions carrying favorable haplotypes and show-
ing a mean score of frost damage ≤1 were mainly winter
fodder peas (e.g. Black seeded, Champagne, Melrose, Blixt
7) (Additional files 6 and 10). While those carrying unfavor-
able haplotypes with a mean score of frost damage ≥4 were
mainly spring garden peas (e.g. Automobile, Caroubel,
Cennia, Ersling) (Additional files 6 and 10). The sequences
of the 75 SNP markers related to frost tolerance trait are
provided in Tayeh et al. [36] (Table S2).

Candidate genes

The projection of the 75 FD-related markers on the pea
genome assembly [39] allowed to define intervals of 2
Mb on all the pea chromosomes, except the chromo-
some 4 (LG IV). Four FD-related markers were assigned
to unanchored scaffolds which were all less than 2 Mb
long: in that case, annotated genes were listed for the
whole scaffold. A particular case was encountered for
the associated marker PsCam036704_21832_970 which
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Quantile-quantile plots of the association mapping results. The plots show the observed p-values (p) for association between a phenotypical trait and each tested marker, expressed as -log10 of p (y-axis) plotted against -log10 of the expected p-values (x-axis) under the null hypothesis of no association for the analyses. **a** is the plot for the trait corresponding to the evaluation of frost damages (FD) in controlled conditions experiment. **b**, **c**, **d** and **e** are the plots for the four traits corresponding to the evaluations of FD in the field experiment</figcaption></figure>

is located on LGVI (51.1 cM) of the genetic consensus map [36] but which projection on the pea genome assembly locates on the chromosome 7 corresponding to the LG VII (Table 3). Considering this ambiguous position, we listed the gene corresponding to this marker as a unique candidate.

We located a total of 867 annotated genes, among which 277 corresponded to genes with unknown function according to the pea genome assembly v.1a (Additional file 11). Among the remaining 590 genes, we focused on gene

families pointed as candidates in previous studies, i.e. *CBF*/*DREB* genes, genes coding for brassinosteroid receptors, genes implied in the production of gibberellin and genes implied in the synthesis of soluble sugars.

Nine candidates corresponding to, or at the vicinity of, FD-related markers were found to be annotated as AP2 domain genes in pea and this annotation could in some cases be refined thanks to the annotation of homologous genes in *M. truncatula* (Additional file 11). The marker PsCam037030_22140_221 on LG VI (49.1 cM), which is
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Quantile-quantile plots of the association mapping results. The plots show the observed p-values (p) for association between a phenotypical trait and each tested marker, expressed as -log10 of p (y-axis) plotted against -log10 of the expected p-values (x-axis) under the null hypothesis of no association for the analyses. **a** is the plot for the trait corresponding to the evaluation of frost damages (FD) in controlled conditions experiment. **b**, **c**, **d** and **e** are the plots for the four traits corresponding to the evaluations of FD in the field experiment</figcaption></figure>

is located on LGVI (51.1 cM) of the genetic consensus map [36] but which projection on the pea genome assembly locates on the chromosome 7 corresponding to the LG VII (Table 3). Considering this ambiguous position, we listed the gene corresponding to this marker as a unique candidate.

We located a total of 867 annotated genes, among which 277 corresponded to genes with unknown function according to the pea genome assembly v.1a (Additional file 11). Among the remaining 590 genes, we focused on gene

families pointed as candidates in previous studies, i.e. *CBF*/*DREB* genes, genes coding for brassinosteroid receptors, genes implied in the production of gibberellin and genes implied in the synthesis of soluble sugars.

Nine candidates corresponding to, or at the vicinity of, FD-related markers were found to be annotated as AP2 domain genes in pea and this annotation could in some cases be refined thanks to the annotation of homologous genes in *M. truncatula* (Additional file 11). The marker PsCam037030_22140_221 on LG VI (49.1 cM), which is
<table><caption>**Table 3** Significant associations detected in the association mapping collection (363 accessions) for the five observed traits</caption><thead><tr><th rowspan="2">Associated marker</th><th colspan="2">Genetic map <sup>a</sup></th><th colspan="2">Physical map <sup>b</sup></th><th rowspan="2">Associated trait (p)</th><th rowspan="2">-log<sub>10</sub> (p)</th><th rowspan="2">Minor allele effect (MAF) <sup>c</sup></th><th rowspan="2">Allelic effect</th><th rowspan="2">SE <sup>e</sup></th><th rowspan="2">Favorable allele</th><th colspan="4">Allele of reference accessions</th></tr><tr><th>Linkage group</th><th>Genetic position (cM)</th><th>Chromosome</th><th>Physical position (bp)</th><th>Champagne</th><th>Térése</th><th>China</th><th>Caméor</th></tr></thead><tbody><tr><td>PsCam020765_11567_1428</td><td>I</td><td>47.5</td><td>2</td><td>160,946,678</td><td>FD Tes Date1</td><td>5.57</td><td>B (0.13)</td><td>0.18</td><td>0.04</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>AA</td></tr><tr><td>PsCam043698_27686_2011</td><td>II</td><td>20.5</td><td>6</td><td>29,326,358</td><td>FD_C_C</td><td>5.45</td><td>B (0.39)</td><td>0.18</td><td>0.04</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>AA</td></tr><tr><td>PsCam049061_31705_537</td><td>III</td><td>131.1</td><td>5</td><td>566,158,491</td><td>FD_C_C</td><td>5.55</td><td>B (0.50)</td><td>-0.29</td><td>0.06</td><td>B</td><td>AA</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam035617_20792_637</td><td>III</td><td>132</td><td>5</td><td>566,953,569</td><td>FD_C_C</td><td>6.63</td><td>A (0.50)</td><td>-0.37</td><td>0.07</td><td>A</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td rowspan="2">PsCam037922_22979_691</td><td rowspan="2">III</td><td rowspan="2">132.1</td><td rowspan="2">5</td><td rowspan="2">567,367,194</td><td>FD Tes Date1</td><td>5.47</td><td>A (0.50)</td><td>-0.31</td><td>0.07</td><td></td><td></td><td></td><td></td></tr><tr><td>FD Tes Date1</td><td>6.63</td><td>A (0.50)</td><td>-0.37</td><td>0.07</td><td>A</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam051352_33858_733</td><td>V</td><td>53.6</td><td>scaffold01940</td><td>189,591</td><td>FD Tes Date4</td><td>7.29</td><td>A (0.25)</td><td>-0.27</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam048068_30823_2326</td><td>V</td><td>56.9</td><td>3</td><td>181,347,340</td><td>FD Tes Date4</td><td>6.21</td><td>B (0.33)</td><td>-0.26</td><td>0.05</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam002670_2123_708</td><td>V</td><td>58.6</td><td>3</td><td>184,294,454</td><td>FD Tes Date4</td><td>6.99</td><td>A (0.19)</td><td>-0.23</td><td>0.04</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam051635_34107_368</td><td>V</td><td>58.6</td><td>3</td><td>184,298,068</td><td>FD Tes Date4</td><td>9.21</td><td>A (0.27)</td><td>-0.33</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam011361_7735_182</td><td>V</td><td>60.2</td><td>3</td><td>189,446,769</td><td>FD Tes Date4</td><td>5.88</td><td>A (0.26)</td><td>-0.25</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam012913_8715_718</td><td>V</td><td>60.2</td><td>3</td><td>187,888,138</td><td>FD Tes Date4</td><td>5.43</td><td>A (0.32)</td><td>-0.22</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam049838_32453_1430</td><td>V</td><td>60.2</td><td>3</td><td>189,010,331</td><td>FD Tes Date4</td><td>7.17</td><td>B (0.28)</td><td>-0.26</td><td>0.05</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam042222_26315_1205</td><td>V</td><td>61</td><td>3</td><td>193,904,545</td><td>FD Tes Date4</td><td>10.01</td><td>A (0.24)</td><td>-0.38</td><td>0.06</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam021402_12023_752</td><td>VI</td><td>41.4</td><td>1</td><td>119,086,683</td><td>FD Tes Date4</td><td>5.84</td><td>A (0.21)</td><td>-0.23</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam043553_27552_265</td><td>VI</td><td>43</td><td>1</td><td>124,617,110</td><td>FD Tes Date4</td><td>5.95</td><td>A (0.22)</td><td>-0.25</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam022275_12527_210</td><td>VI</td><td>47.4</td><td>1</td><td>167,098,866</td><td>FD_C_C</td><td>6.21</td><td>B (0.23)</td><td>-0.33</td><td>0.06</td><td>B</td><td>BB</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td rowspan="2">PsCam003541_2731_574</td><td rowspan="2">VI</td><td rowspan="2">48</td><td rowspan="2">1</td><td rowspan="2">165,049,951</td><td>FD Tes Date1</td><td>5.64</td><td>B (0.23)</td><td>-0.29</td><td>0.05</td><td></td><td></td><td></td><td></td></tr><tr><td>FD Tes Date4</td><td>6.86</td><td>B (0.23)</td><td>-0.27</td><td>0.05</td><td></td><td></td><td></td><td></td></tr><tr><td>PsCam050834_33394_133</td><td>VI</td><td>48.5</td><td>1</td><td>165,743,845</td><td>FD_C_C</td><td>6.51</td><td>A (0.44)</td><td>0.21</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>BB</td></tr><tr><td>PsCam001662_1375_798</td><td>VI</td><td>49.1</td><td>1</td><td>168,461,782</td><td>FD_C_C</td><td>5.35</td><td>A (0.33)</td><td>-0.25</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam043847_27818_181</td><td>VI</td><td>49.1</td><td>1</td><td>175,268,097</td><td>FD Tes Date4</td><td>6.08</td><td>A (0.31)</td><td>-0.24</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam045419_29091_815</td><td>VI</td><td>49.1</td><td>1</td><td>169,979,594</td><td>FD_C_C</td><td>7.57</td><td>B (0.37)</td><td>-0.24</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam048404_31138_2325</td><td>VI</td><td>49.1</td><td>1</td><td>170,956,086</td><td>FD_C_C</td><td>6.02</td><td>B (0.31)</td><td>-0.30</td><td>0.06</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam004890_3696_1580</td><td>VI</td><td>49.1</td><td>1</td><td>181,902,508</td><td>FD Tes Date4</td><td>5.82</td><td>A (0.22)</td><td>-0.24</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam023246_13111_1125</td><td>VI</td><td>49.1</td><td>1</td><td>167,654,791</td><td>FD Tes Date4</td><td>5.40</td><td>B (0.21)</td><td>-0.25</td><td>0.05</td><td>B</td><td>BB</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td>PsCam037082_22189_1302</td><td>VI</td><td>49.1</td><td>1</td><td>181,965,322</td><td>FD Tes Date4</td><td>5.40</td><td>B (0.22)</td><td>-0.23</td><td>0.05</td><td>B</td><td>BB</td><td>AA</td><td>AA</td><td>AA</td></tr></tbody></table>
<table><caption>**Table 3** Significant associations detected in the association mapping collection (363 accessions) for the five observed traits</caption><thead><tr><th rowspan="2">Associated marker</th><th colspan="2">Genetic map <sup>a</sup></th><th colspan="2">Physical map <sup>b</sup></th><th rowspan="2">Associated trait (p)</th><th rowspan="2">-log<sub>10</sub> (p)</th><th rowspan="2">Minor allele effect (MAF) <sup>c</sup></th><th rowspan="2">Allelic effect</th><th rowspan="2">SE <sup>e</sup></th><th rowspan="2">Favorable allele</th><th colspan="3">Allele of reference accessions</th></tr><tr><th>Linkage group</th><th>Genetic position (cM)</th><th>Chromosome</th><th>Physical position (bp)</th><th>Champagne</th><th>Térése</th><th>China</th><th>Caméor</th></tr></thead><tbody><tr><td>PsCam020765_11567_1428</td><td>I</td><td>47.5</td><td>2</td><td>160,946,678</td><td>FD Tes Date1</td><td>5.57</td><td>B (0.13)</td><td>0.18</td><td>0.04</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>AA</td></tr><tr><td>PsCam043698_27686_2011</td><td>II</td><td>20.5</td><td>6</td><td>29,326,358</td><td>FD_C_C</td><td>5.45</td><td>B (0.39)</td><td>0.18</td><td>0.04</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>AA</td></tr><tr><td>PsCam049061_31705_537</td><td>III</td><td>131.1</td><td>5</td><td>566,158,491</td><td>FD_C_C</td><td>5.55</td><td>B (0.50)</td><td>-0.29</td><td>0.06</td><td>B</td><td>AA</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam035617_20792_637</td><td>III</td><td>132</td><td>5</td><td>566,953,569</td><td>FD_C_C</td><td>6.63</td><td>A (0.50)</td><td>-0.37</td><td>0.07</td><td>A</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td rowspan="2">PsCam037922_22979_691</td><td rowspan="2">III</td><td rowspan="2">132.1</td><td rowspan="2">5</td><td rowspan="2">567,367,194</td><td>FD Tes Date1</td><td>5.47</td><td>A (0.50)</td><td>-0.31</td><td>0.07</td><td></td><td></td><td></td><td></td></tr><tr><td>FD Tes Date1</td><td>6.63</td><td>A (0.50)</td><td>-0.37</td><td>0.07</td><td>A</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam051352_33858_733</td><td>V</td><td>53.6</td><td>scaffold01940</td><td>189,591</td><td>FD Tes Date4</td><td>7.29</td><td>A (0.25)</td><td>-0.27</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam048068_30823_2326</td><td>V</td><td>56.9</td><td>3</td><td>181,347,340</td><td>FD Tes Date4</td><td>6.21</td><td>B (0.33)</td><td>-0.26</td><td>0.05</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam002670_2123_708</td><td>V</td><td>58.6</td><td>3</td><td>184,294,454</td><td>FD Tes Date4</td><td>6.99</td><td>A (0.19)</td><td>-0.23</td><td>0.04</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam051635_34107_368</td><td>V</td><td>58.6</td><td>3</td><td>184,298,068</td><td>FD Tes Date4</td><td>9.21</td><td>A (0.27)</td><td>-0.33</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam011361_7735_182</td><td>V</td><td>60.2</td><td>3</td><td>189,446,769</td><td>FD Tes Date4</td><td>5.88</td><td>A (0.26)</td><td>-0.25</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam012913_8715_718</td><td>V</td><td>60.2</td><td>3</td><td>187,888,138</td><td>FD Tes Date4</td><td>5.43</td><td>A (0.32)</td><td>-0.22</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam049838_32453_1430</td><td>V</td><td>60.2</td><td>3</td><td>189,010,331</td><td>FD Tes Date4</td><td>7.17</td><td>B (0.28)</td><td>-0.26</td><td>0.05</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam042222_26315_1205</td><td>V</td><td>61</td><td>3</td><td>193,904,545</td><td>FD Tes Date4</td><td>10.01</td><td>A (0.24)</td><td>-0.38</td><td>0.06</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam021402_12023_752</td><td>VI</td><td>41.4</td><td>1</td><td>119,086,683</td><td>FD Tes Date4</td><td>5.84</td><td>A (0.21)</td><td>-0.23</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam043553_27552_265</td><td>VI</td><td>43</td><td>1</td><td>124,617,110</td><td>FD Tes Date4</td><td>5.95</td><td>A (0.22)</td><td>-0.25</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam022275_12527_210</td><td>VI</td><td>47.4</td><td>1</td><td>167,098,866</td><td>FD_C_C</td><td>6.21</td><td>B (0.23)</td><td>-0.33</td><td>0.06</td><td>B</td><td>BB</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td rowspan="2">PsCam003541_2731_574</td><td rowspan="2">VI</td><td rowspan="2">48</td><td rowspan="2">1</td><td rowspan="2">165,049,951</td><td>FD Tes Date1</td><td>5.64</td><td>B (0.23)</td><td>-0.29</td><td>0.05</td><td></td><td></td><td></td><td></td></tr><tr><td>FD Tes Date4</td><td>6.86</td><td>B (0.23)</td><td>-0.27</td><td>0.05</td><td></td><td></td><td></td><td></td></tr><tr><td>PsCam050834_33394_133</td><td>VI</td><td>48.5</td><td>1</td><td>165,743,845</td><td>FD_C_C</td><td>6.51</td><td>A (0.44)</td><td>0.21</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>BB</td></tr><tr><td>PsCam001662_1375_798</td><td>VI</td><td>49.1</td><td>1</td><td>168,461,782</td><td>FD_C_C</td><td>5.35</td><td>A (0.33)</td><td>-0.25</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam043847_27818_181</td><td>VI</td><td>49.1</td><td>1</td><td>175,268,097</td><td>FD Tes Date4</td><td>6.08</td><td>A (0.31)</td><td>-0.24</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam045419_29091_815</td><td>VI</td><td>49.1</td><td>1</td><td>169,979,594</td><td>FD_C_C</td><td>5.90</td><td>A (0.27)</td><td>-0.29</td><td>0.06</td><td>A</td><td>AA</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam048404_31138_2325</td><td>VI</td><td>49.1</td><td>1</td><td>170,956,086</td><td>FD_C_C</td><td>7.57</td><td>B (0.37)</td><td>-0.24</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam004890_3696_1580</td><td>VI</td><td>49.1</td><td>1</td><td>181,902,508</td><td>FD Tes Date4</td><td>6.02</td><td>B (0.31)</td><td>-0.30</td><td>0.06</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam023246_13111_1125</td><td>VI</td><td>49.1</td><td>1</td><td>167,654,791</td><td>FD Tes Date4</td><td>5.82</td><td>A (0.22)</td><td>-0.24</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam037082_22189_1302</td><td>VI</td><td>49.1</td><td>1</td><td>181,965,322</td><td>FD Tes Date4</td><td>5.40</td><td>B (0.22)</td><td>-0.23</td><td>0.05</td><td>B</td><td>BB</td><td>AA</td><td>AA</td><td>AA</td></tr></tbody></table>
<table><caption>**Table 3** Significant associations detected in the association mapping collection (363 accessions) for the five observed traits (Continued)</caption><thead><tr><th rowspan="2">Associated marker</th><th colspan="2">Genetic map <sup>a</sup></th><th colspan="2">Physical map <sup>b</sup></th><th rowspan="2">Associated trait (p)</th><th colspan="2">-log<sub>10</sub> (p)</th><th rowspan="2">Minor allele (MAF) <sup>c</sup></th><th rowspan="2">Allelic effect <sup>d</sup></th><th rowspan="2">SE <sup>e</sup></th><th rowspan="2">Favorable allele</th><th colspan="3">Allele of reference accessions</th></tr><tr><th>Linkage group</th><th>Genetic position (cM)</th><th>Chromosome</th><th>Physical position (bp)</th><th>Champagne</th><th>Térése</th><th>China</th><th>Caméor</th></tr></thead><tbody><tr><td>PsCam053556_35439_228</td><td>VI</td><td>49.1</td><td>1</td><td>166,361,064</td><td>FD Tes date4</td><td>5.39</td><td>B (0.31)</td><td>-0.21</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam007060_5248_2156</td><td>VI</td><td>50.1</td><td>1</td><td>179,745,943</td><td>FD_Cc</td><td>7.66</td><td>B (0.40)</td><td>0.22</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td><td>BB</td></tr><tr><td>PsCam026381_15234_450</td><td>VI</td><td>50.1</td><td>1</td><td>178,695,180</td><td>FD_Cc</td><td>6.33</td><td>B (0.47)</td><td>0.21</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>BB</td></tr><tr><td>PsCam037879_22937_948</td><td>VI</td><td>50.1</td><td>1</td><td>178,713,424</td><td>FD_Cc</td><td>9.90</td><td>A (0.43)</td><td>0.26</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam001302_1089_1249</td><td>VI</td><td>50.2</td><td>1</td><td>182,819,711</td><td>FD_Cc</td><td>6.50</td><td>B (0.46)</td><td>-0.22</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>AA</td></tr><tr><td>PsCam003337_2589_220</td><td>VI</td><td>50.2</td><td>1</td><td>180,428,730</td><td>FD_Cc</td><td>7.34</td><td>A (0.39)</td><td>0.22</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam037603_22679_2024</td><td>VI</td><td>50.2</td><td>1</td><td>180,042,543</td><td>FD_Cc</td><td>7.59</td><td>A (0.38)</td><td>0.22</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam040539_25267_719</td><td>VI</td><td>50.2</td><td>1</td><td>181,440,108</td><td>FD_Cc</td><td>7.59</td><td>A (0.38)</td><td>0.22</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>AA</td></tr><tr><td>PsCam035398_20586_1094</td><td>VI</td><td>50.6</td><td>1</td><td>194,261,395</td><td>FD_Cc</td><td>7.63</td><td>B (0.50)</td><td>-0.25</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam039054_24037_586</td><td>VI</td><td>50.6</td><td>1</td><td>210,266,228</td><td>FD_Cc</td><td>6.30</td><td>B (0.41)</td><td>0.21</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td><td>BB</td></tr><tr><td>PsCam049925_32537_2703</td><td>VI</td><td>50.7</td><td>1</td><td>194,874,039</td><td>FD_Cc</td><td>6.57</td><td>B (0.41)</td><td>0.22</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td><td>BB</td></tr><tr><td>PsCam038477_23510_469</td><td>VI</td><td>50.8</td><td>1</td><td>194,779,560</td><td>FD_Cc</td><td>7.67</td><td>B (0.39)</td><td>-0.27</td><td>0.05</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam001148_968_692</td><td>VI</td><td>51.1</td><td>1</td><td>218,959,228</td><td>FD_Cc</td><td>6.06</td><td>B (0.50)</td><td>0.22</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td>PsCam006586_4888_245</td><td>VI</td><td>51.1</td><td>1</td><td>218,476,070</td><td>FD_Cc</td><td>6.00</td><td>A (0.41)</td><td>0.20</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam017327_10745_93</td><td>VI</td><td>51.1</td><td>1</td><td>208,681,226</td><td>FD_Cc</td><td>8.40</td><td>B (0.48)</td><td>-0.27</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>AA</td></tr><tr><td>PsCam020872_11644_1756</td><td>VI</td><td>51.1</td><td>1</td><td>217,751,686</td><td>FD_Cc</td><td>7.31</td><td>B (0.45)</td><td>0.23</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td>PsCam025448_14564_2314</td><td>VI</td><td>51.1</td><td>1</td><td>218,939,803</td><td>FD_Cc</td><td>6.08</td><td>A (0.29)</td><td>-0.22</td><td>0.04</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam033562_19160_307</td><td>VI</td><td>51.1</td><td>1</td><td>217,753,744</td><td>FD_Cc</td><td>6.04</td><td>A (0.48)</td><td>0.21</td><td>0.04</td><td>B</td><td>AA</td><td>AA</td><td>BB</td><td>BB</td></tr><tr><td>PsCam036704_21832_970</td><td>VI</td><td>51.1</td><td>7</td><td>404,463,902</td><td>FD_Cc</td><td>7.12</td><td>B (0.48)</td><td>0.22</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam038677_23691_325</td><td>VI</td><td>51.1</td><td>1</td><td>221,377,153</td><td>FD_Cc</td><td>6.59</td><td>B (0.43)</td><td>0.20</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>BB</td></tr><tr><td>PsCam044418_28295_579</td><td>VI</td><td>51.1</td><td>1</td><td>217,447,863</td><td>FD_Cc</td><td>6.61</td><td>A (0.41)</td><td>0.21</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam052151_34533_410</td><td>VI</td><td>51.1</td><td>1</td><td>216,779,912</td><td>FD_Cc</td><td>6.16</td><td>A (0.47)</td><td>0.21</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>BB</td></tr><tr><td>PsCam053189_35238_569</td><td>VI</td><td>51.1</td><td>1</td><td>216,188,812</td><td>FD_Cc</td><td>7.30</td><td>B (0.49)</td><td>0.23</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam042568_26621_42</td><td>VI</td><td>51.1</td><td>1</td><td>213,838,599</td><td>FD Tes date4</td><td>5.43</td><td>A (0.21)</td><td>-0.25</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam035356_20546_778</td><td>VI</td><td>52.5</td><td>1</td><td>221,502,038</td><td>FD_Cc</td><td>5.99</td><td>B (0.43)</td><td>0.20</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam040326_25099_1923</td><td>VI</td><td>52.5</td><td>1</td><td>226,457,822</td><td>FD_Cc</td><td>5.91</td><td>A (0.38)</td><td>0.20</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam011774_8038_200</td><td>VI</td><td>52.8</td><td>1</td><td>232,110,795</td><td>FD_Cc</td><td>6.76</td><td>A (0.49)</td><td>-0.23</td><td>0.04</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam023852_13541_264</td><td>VI</td><td>52.8</td><td>1</td><td>232,209,522</td><td>FD_Cc</td><td>6.71</td><td>B (0.46)</td><td>-0.24</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>AA</td><td>BB</td></tr><tr><td>PsCam026055_15000_192</td><td>VI</td><td>52.8</td><td>1</td><td>228,344,153</td><td>FD_Cc</td><td>5.70</td><td>B (0.44)</td><td>-0.22</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>AA</td><td>BB</td></tr><tr><td>PsCam001108_940_48</td><td>VII</td><td>72.9</td><td>7</td><td>338,172,762</td><td>FD Tes date2</td><td>5.74</td><td>B (0.48)</td><td>0.20</td><td>0.04</td><td>A</td><td>BB</td><td>AA</td><td>AA</td><td>AA</td></tr></tbody></table>
<table><caption>**Table 3** Significant associations detected in the association mapping collection (363 accessions) for the five observed traits (Continued)</caption><thead><tr><th rowspan="2">Associated marker</th><th colspan="2">Genetic map <sup>a</sup></th><th colspan="2">Physical map <sup>b</sup></th><th rowspan="2">Associated trait (p)</th><th colspan="2">-log<sub>10</sub> (p)</th><th colspan="2">Allelic effect</th><th rowspan="2">SE <sup>c</sup></th><th rowspan="2">Favorable allele</th><th colspan="2">Allele of reference accessions</th></tr><tr><th>Linkage group</th><th>Genetic position (cM)</th><th>Chromosome</th><th>Physical position (bp)</th><th>Minor allele (MAF)</th><th>Allelic effect</th><th>Champagne</th><th>Térése</th><th>China</th><th>Caméor</th></tr></thead><tbody><tr><td>PsCam053556_35439_228</td><td>VI</td><td>49.1</td><td>1</td><td>166,361,064</td><td>FD Tes date4</td><td>5.39</td><td>B (0.31)</td><td>-0.21</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam007060_5248_2156</td><td>VI</td><td>50.1</td><td>1</td><td>179,745,943</td><td>FD_C_C</td><td>7.66</td><td>B (0.40)</td><td>0.22</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td><td>BB</td></tr><tr><td>PsCam026381_15234_450</td><td>VI</td><td>50.1</td><td>1</td><td>178,695,180</td><td>FD_C_C</td><td>6.33</td><td>B (0.47)</td><td>0.21</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>BB</td></tr><tr><td>PsCam037879_22937_948</td><td>VI</td><td>50.1</td><td>1</td><td>178,713,424</td><td>FD_C_C</td><td>9.90</td><td>A (0.43)</td><td>0.26</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam001302_1089_1249</td><td>VI</td><td>50.2</td><td>1</td><td>182,819,711</td><td>FD_C_C</td><td>6.50</td><td>B (0.46)</td><td>-0.22</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>AA</td></tr><tr><td>PsCam003337_2589_220</td><td>VI</td><td>50.2</td><td>1</td><td>180,428,730</td><td>FD_C_C</td><td>7.34</td><td>A (0.39)</td><td>0.22</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam037603_22679_2024</td><td>VI</td><td>50.2</td><td>1</td><td>180,042,543</td><td>FD_C_C</td><td>7.59</td><td>A (0.38)</td><td>0.22</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam040539_25267_719</td><td>VI</td><td>50.2</td><td>1</td><td>181,440,108</td><td>FD_C_C</td><td>7.59</td><td>A (0.38)</td><td>0.22</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>AA</td></tr><tr><td>PsCam035398_20586_1094</td><td>VI</td><td>50.6</td><td>1</td><td>194,261,395</td><td>FD_C_C</td><td>7.63</td><td>B (0.50)</td><td>-0.25</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam039054_24037_586</td><td>VI</td><td>50.6</td><td>1</td><td>210,266,228</td><td>FD_C_C</td><td>6.30</td><td>B (0.41)</td><td>0.21</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td><td>BB</td></tr><tr><td>PsCam049925_32537_2703</td><td>VI</td><td>50.7</td><td>1</td><td>194,874,039</td><td>FD_C_C</td><td>6.57</td><td>B (0.41)</td><td>0.22</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td><td>BB</td></tr><tr><td>PsCam038477_23510_469</td><td>VI</td><td>50.8</td><td>1</td><td>194,779,560</td><td>FD_C_C</td><td>7.67</td><td>B (0.39)</td><td>-0.27</td><td>0.05</td><td>B</td><td>BB</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam001148_968_692</td><td>VI</td><td>51.1</td><td>1</td><td>218,959,228</td><td>FD_C_C</td><td>6.06</td><td>B (0.50)</td><td>0.22</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam006586_4888_245</td><td>VI</td><td>51.1</td><td>1</td><td>218,476,070</td><td>FD_C_C</td><td>6.00</td><td>A (0.41)</td><td>0.20</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam017327_10745_93</td><td>VI</td><td>51.1</td><td>1</td><td>208,681,226</td><td>FD_C_C</td><td>8.40</td><td>B (0.48)</td><td>-0.27</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>AA</td></tr><tr><td>PsCam020872_11644_1756</td><td>VI</td><td>51.1</td><td>1</td><td>217,751,686</td><td>FD_C_C</td><td>7.31</td><td>B (0.45)</td><td>0.23</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam025448_14564_2314</td><td>VI</td><td>51.1</td><td>1</td><td>218,939,803</td><td>FD_C_C</td><td>6.08</td><td>A (0.29)</td><td>-0.22</td><td>0.04</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam033562_19160_307</td><td>VI</td><td>51.1</td><td>1</td><td>217,753,744</td><td>FD_C_C</td><td>6.04</td><td>A (0.48)</td><td>0.21</td><td>0.04</td><td>B</td><td>AA</td><td>AA</td><td>BB</td><td>BB</td></tr><tr><td>PsCam036704_21832_970</td><td>VI</td><td>51.1</td><td>7</td><td>404,463,902</td><td>FD_C_C</td><td>7.12</td><td>B (0.48)</td><td>0.22</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam038677_23691_325</td><td>VI</td><td>51.1</td><td>1</td><td>221,377,153</td><td>FD_C_C</td><td>6.59</td><td>B (0.43)</td><td>0.20</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>BB</td></tr><tr><td>PsCam044418_28295_579</td><td>VI</td><td>51.1</td><td>1</td><td>217,447,863</td><td>FD_C_C</td><td>6.61</td><td>A (0.41)</td><td>0.21</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam052151_34533_410</td><td>VI</td><td>51.1</td><td>1</td><td>216,779,912</td><td>FD_C_C</td><td>6.16</td><td>A (0.47)</td><td>0.21</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>BB</td></tr><tr><td>PsCam053189_35238_569</td><td>VI</td><td>51.1</td><td>1</td><td>216,188,812</td><td>FD_C_C</td><td>7.30</td><td>B (0.49)</td><td>0.23</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam042568_26621_42</td><td>VI</td><td>51.1</td><td>1</td><td>213,838,599</td><td>FD Tes date4</td><td>5.43</td><td>A (0.21)</td><td>-0.25</td><td>0.05</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>BB</td></tr><tr><td>PsCam035356_20546_778</td><td>VI</td><td>52.5</td><td>1</td><td>221,502,038</td><td>FD_C_C</td><td>5.99</td><td>B (0.43)</td><td>0.20</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>BB</td><td>AA</td></tr><tr><td>PsCam040326_25099_1923</td><td>VI</td><td>52.5</td><td>1</td><td>226,457,822</td><td>FD_C_C</td><td>5.91</td><td>A (0.38)</td><td>0.20</td><td>0.04</td><td>B</td><td>BB</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam011774_8038_200</td><td>VI</td><td>52.8</td><td>1</td><td>232,110,795</td><td>FD_C_C</td><td>6.76</td><td>A (0.49)</td><td>-0.23</td><td>0.04</td><td>A</td><td>AA</td><td>BB</td><td>BB</td><td>AA</td></tr><tr><td>PsCam023852_13541_264</td><td>VI</td><td>52.8</td><td>1</td><td>232,209,522</td><td>FD_C_C</td><td>6.71</td><td>B (0.46)</td><td>-0.24</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>AA</td><td>BB</td></tr><tr><td>PsCam026055_15000_192</td><td>VI</td><td>52.8</td><td>1</td><td>228,344,153</td><td>FD_C_C</td><td>5.70</td><td>B (0.44)</td><td>-0.22</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>AA</td><td>BB</td></tr><tr><td>PsCam001108_940_48</td><td>VII</td><td>72.9</td><td>7</td><td>338,172,762</td><td>FD Tes date2</td><td>5.74</td><td>B (0.48)</td><td>0.20</td><td>0.04</td><td>A</td><td>BB</td><td>AA</td><td>AA</td><td>AA</td></tr></tbody></table>
<table><caption>**Table 3** Significant associations detected in the association mapping collection (363 accessions) for the five observed traits (Continued)</caption><thead><tr><th rowspan="2">Associated marker</th><th colspan="2">Genetic map <sup>a</sup></th><th colspan="2">Physical map <sup>b</sup></th><th rowspan="2">Associated trait (p)</th><th rowspan="2">-log<sub>10</sub> (p)</th><th rowspan="2">Minor allele (MAF) <sup>c</sup></th><th rowspan="2">Allelic effect <sup>d</sup></th><th rowspan="2">SE <sup>e</sup></th><th rowspan="2">Favorable allele</th><th colspan="3">Allele of reference accessions</th></tr><tr><th>Linkage group</th><th>Genetic position (cM)</th><th>Chromosome</th><th>Physical position (bp)</th><th>Champagne</th><th>Térèse</th><th>China</th><th>Caméor</th></tr></thead><tbody><tr><td>PsCam050827_33387_242</td><td>VII</td><td>82.4</td><td>scaffold02790</td><td>184,253</td><td>FD Tes Date3</td><td>5.34</td><td>B (0.48)</td><td>0.21</td><td>0.05</td><td></td><td></td><td></td><td></td></tr><tr><td>PsCam037927_22984_97</td><td>VII</td><td>82.5</td><td>7</td><td>364,561,668</td><td>FD Tes Date2</td><td>6.27</td><td>B (0.21)</td><td>0.21</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td>PsCam038378_23415_721</td><td>VII</td><td>86.4</td><td>7</td><td>416,077,744</td><td>FD CC</td><td>5.40</td><td>B (0.41)</td><td>0.19</td><td>0.04</td><td>A</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam042427_26489_524</td><td>VII</td><td>87.1</td><td>7</td><td>417,288,899</td><td>FD CC</td><td>6.18</td><td>A (0.46)</td><td>0.22</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam000301_265_1519</td><td>VII</td><td>87.9</td><td>7</td><td>424,418,432</td><td>FD Tes Date2</td><td>5.67</td><td>B (0.41)</td><td>0.20</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td>PsCam045010_28743_3365</td><td>VII</td><td>88</td><td>scaffold02266</td><td>42,102</td><td>FD Tes Date3</td><td>7.41</td><td>B (0.19)</td><td>0.23</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td>PsCam004928_3732_3087</td><td>VII</td><td>89.3</td><td>7</td><td>430,325,609</td><td>FD Tes Date2</td><td>5.72</td><td>B (0.19)</td><td>0.22</td><td>0.05</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>FD Tes Date2</td><td>5.73</td><td>B (0.21)</td><td>0.20</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>FD Tes Date2</td><td>7.46</td><td>B (0.28)</td><td>0.25</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>FD Tes Date3</td><td>5.40</td><td>B (0.28)</td><td>0.24</td><td>0.05</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Traits are coded as described in Table 1. For each associated SNP, genetic and physical positions are indicated. <sup>a</sup>: Linkage groups are named from I to VII and the genetic position is indicated in cM Haldane following the consensus map from Tayeh et al. [36]. <sup>b</sup>: Physical marker position results from the projection of the marker onto the *Pisum sativum* v.1a genome JBrowse available at https://urgi.versailles.inra.fr/Species/Pisum [39]. Chromosomes are named in consecutive numerical order from 1 to 7 and physical positions are indicated in bp. Markers which are not assigned to one of the seven chromosomes are represented by their physical position on unanchored scaffolds. Significance threshold of the marker-trait association is -log₁₀ (p) > 5.33 as described in the Material and Methods section. For each associated SNP, the allelic effect of the minor allele as well as the favorable allele are listed. <sup>c</sup>: Minor allele frequency; <sup>d</sup>: Allele effect of the minor allele. <sup>e</sup>: standard error of the allelic effect. Alleles of four cultivars which are the progenitors of bi-parental mapping populations for frost tolerance are listed to illustrate contrasted genotypes for this trait: Champagne and China (frost tolerant), Térèse and Caméor (frost sensitive)
<table><caption>**Table 3** Significant associations detected in the association mapping collection (363 accessions) for the five observed traits (Continued)</caption><thead><tr><th rowspan="2">Associated marker</th><th colspan="2">Genetic map <sup>a</sup></th><th colspan="2">Physical map <sup>b</sup></th><th rowspan="2">Associated trait (p)</th><th rowspan="2">-log<sub>10</sub> (p)</th><th rowspan="2">Minor allele (MAF) <sup>c</sup></th><th rowspan="2">Allelic effect <sup>d</sup></th><th rowspan="2">SE <sup>e</sup></th><th rowspan="2">Favorable allele</th><th colspan="3">Allele of reference accessions</th></tr><tr><th>Linkage group</th><th>Genetic position (cM)</th><th>Chromosome</th><th>Physical position (bp)</th><th>Champagne</th><th>Térèse</th><th>China</th><th>Caméor</th></tr></thead><tbody><tr><td>PsCam050827_33387_242</td><td>VII</td><td>82.4</td><td>scaffold02790</td><td>184,253</td><td>FD Tes Date3</td><td>5.34</td><td>B (0.48)</td><td>0.21</td><td>0.05</td><td></td><td></td><td></td><td></td></tr><tr><td>PsCam037927_22984_97</td><td>VII</td><td>82.5</td><td>7</td><td>364,561,668</td><td>FD Tes Date2</td><td>6.27</td><td>B (0.21)</td><td>0.21</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td>PsCam038378_23415_721</td><td>VII</td><td>86.4</td><td>7</td><td>416,077,744</td><td>FD CC</td><td>5.40</td><td>B (0.41)</td><td>0.19</td><td>0.04</td><td>A</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam042427_26489_524</td><td>VII</td><td>87.1</td><td>7</td><td>417,288,899</td><td>FD CC</td><td>6.18</td><td>A (0.46)</td><td>0.22</td><td>0.04</td><td>B</td><td>BB</td><td>AA</td><td>BB</td></tr><tr><td>PsCam000301_265_1519</td><td>VII</td><td>87.9</td><td>7</td><td>424,418,432</td><td>FD Tes Date2</td><td>5.67</td><td>B (0.41)</td><td>0.20</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td>PsCam045010_28743_3365</td><td>VII</td><td>88</td><td>scaffold02266</td><td>42,102</td><td>FD Tes Date3</td><td>5.72</td><td>B (0.19)</td><td>0.23</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td>PsCam004928_3732_3087</td><td>VII</td><td>89.3</td><td>7</td><td>430,325,609</td><td>FD Tes Date2</td><td>5.73</td><td>B (0.21)</td><td>0.20</td><td>0.04</td><td>A</td><td>AA</td><td>AA</td><td>AA</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>FD Tes Date3</td><td>5.40</td><td>B (0.28)</td><td>0.24</td><td>0.05</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Traits are coded as described in Table 1. For each associated SNP, genetic and physical positions are indicated. <sup>a</sup>: Linkage groups are named from I to VII and the genetic position is indicated in cM Haldane following the consensus map from Tayeh et al. [36]. <sup>b</sup>: Physical marker position results from the projection of the marker onto the *Pisum sativum* v.1a genome JBrowse available at https://urgi.versailles.inra.fr/Species/Pisum [39]. Chromosomes are named in consecutive numerical order from 1 to 7 and physical positions are indicated in bp. Markers which are not assigned to one of the seven chromosomes are represented by their physical position on unanchored scaffolds. Significance threshold of the marker-trait association is -log<sub>10</sub> (p) > 5.33 as described in the Material and Methods section. For each associated SNP, the allelic effect of the minor allele as well as the favorable allele are listed. <sup>c</sup>: Minor allele frequency; <sup>d</sup>: Allele effect of the minor allele. <sup>e</sup>: standard error of the allelic effect. Alleles of four cultivars which are the progenitors of bi-parental mapping populations for frost tolerance are listed to illustrate contrasted genotypes for this trait: Champagne and China (frost tolerant), Térèse and Caméor (frost sensitive)
in high LD with all the significant FD-associated markers belonging to the LD block VI.2, belongs to the gene Psat1g103560 annotated as a *CBF* gene according to the annotation of the homologous gene of *M. truncatula* (Additional file 11). This FD-related marker was also annotated as a *CBF14* gene, following a blast search against *M. truncatula* carried out by Tayeh et al. [36]. Two other genes (Psat1g103560 and Psat1g103600) annotated as *CBF* genes were identified close to PsCam037030_22140_221, one of which being also precisely annotated as a *CBF14* by Tayeh et al. [36]. In the LD block VI.2, 2 FD-associated markers, namely PsCam023246_13111_1125 and PsCam007060_5248_2156, were also found to be close to AP2/ERF (Ethylene-responsive transcription factor) genes, namely Psat1g097280 and Psat1g097280, for which a precise study of the sequence is needed to verify if they belong to the CBF sub-family. Finally, 4 other potential candidate genes were found at the vicinity of the LD block VI.2, namely Psat1g103920, Psat1g106640, Psat1g115640 and Psat1g103680. These genes were annotated as *DREB* (Dehydration-responsive element-binding protein) genes, another term used to refer to *CBF* genes, in the homologous genes of *M. truncatula* (Additional file 11). These 4 genes lie in an interval of 35 Mb situated at a distance of 448 kb from the three other *CBF* genes mentioned above. One of them, Psat1g103920, contained the marker PsCam050192_32788_145 which was excluded from the GWA analysis because it showed a minor allele frequency lower than 0.05.

Three FD-associated markers, namely PsCam035617_20792_637 (LD-block III.1), PsCam048068_30823_2326 (LD-block V.1) and PsCam011774_8038_200 (LD-block VI.8) were found to correspond to 3 genes, namely Psat5g299600, Psat3g087400 and Psat1g119400 encoding brassinosteroid receptors (Additional file 11).

The FD-related marker PsCam037922_22979_691, which is only 0.1 cM apart from PsCam035617_20792_637 just mentioned above, was found to lie in Psat5g299720, a gene encoding for the gibberellin 3beta-hydroxylase enzyme and shown to correspond to the dwarfism gene *Le* in pea (Additional file 11).

Finally, three FD-associated markers belonging to the FD-related locus on LG VII, namely PsCam001108_940_48, PsCam037927_22984_97 and PsCam004928_3732_3087 were located within genes related to the synthesis of soluble sugars (Psat7g180280: endo-beta-1,3-glucanase, Psat7g193120: endo-1,4-beta-glucanase and Psat7g214880: beta-glucosidase G2, respectively) (Additional file 11).

# Discussion

## GWAS brings new insights into the determinism of frost tolerance in pea

In the present study, 75 markers associated with frost tolerance, 62 markers significantly detected by GWAS

and 13 markers in high LD ($r^2 > 0.8$) with one or the other of the 62 markers, were located among all the linkage groups of pea, except LG IV (Table 3, Additional file 9).

Comparing the map positions with those of frost tolerance QTLs previously described, 3 regions corresponding to 62 markers among the 75 markers, were found to colocalize with 3 main QTLs previously detected by linkage mapping in two bi-parental populations for the same trait, namely WFD 3.2, WFD 5.1 and WFD 6.1 [21–23] (Fig. 4). These 3 QTLs were repeatedly detected in 5, 11 and 10 field conditions for WFD 3.2, WFD 5.1 and WFD 6.1, respectively. Moreover, the position corresponding to WFD 6.1 seems to also match with an EST marker recently found to be associated with frost tolerance in pea [38]: indeed, Liu et al. identified 7 marker-trait associations within a collection of 672 accessions, among which the marker EST1109 was located on LG VI within a functional gene that has a high homology with a gene encoding an alpha-mannosidase in *M. truncatula*. We reviewed the 1646 transcript-derived SNP markers mapped on LG VI by Tayeh et al. [36] and found that a marker corresponding to an alpha-mannosidase-like protein was located at 43.2 cM on the consensus map; consequently, as the FD-associated markers detected on LG VI in the present study are comprised between 41.4 and 52.8 cM, it is likely that the LG VI positions of both association studies coincide. Similarly, this study validated the QTL previously detected on LG III (WFD 3.2 in Pop2 [21], III.2 in Pop9 [23]) with a higher resolution than previous linkage mapping studies by the identifications of three main significant FD-associated markers, and whose positive alleles decreased the frost damage of accessions by an average of 0.33 (Table 3). All favorable alleles were carried by the sensitive genotypes Térèse and Caméor unlike the tolerant genotypes Champagne and China which had all undesirable alleles underlying this locus (Table 3). Altogether, the consistency of the 3 positions detected in both bi-parental mapping and association genetics reinforces their interest for breeding. The results presented here constitute an additional step towards the identification of underlying genes potentially involved in the control of frost tolerance, thanks to refined intervals provided by GWAS.

Unlike the correspondences with bi-parental mapping positions presented above, the present GWA study did not highlight any colocalization with a major QTL, namely WFD3.1, which was however found to be responsible for up to 52 and 19% of the winter frost damage variation within the RILs populations derived from Champagne x Térèse (Pop2) and China x Caméor (Pop9), respectively [21, 23]. The flowering gene *Hr* (*High response to photoperiod*), an orthologue of the Arabidopsis *Early flowering 3*, *Elf3* [40], was shown to be
in high LD with all the significant FD-associated markers belonging to the LD block VI.2, belongs to the gene Psat1g103560 annotated as a *CBF* gene according to the annotation of the homologous gene of *M. truncatula* (Additional file 11). This FD-related marker was also annotated as a *CBF14* gene, following a blast search against *M. truncatula* carried out by Tayeh et al. [36]. Two other genes (Psat1g103560 and Psat1g103600) annotated as *CBF* genes were identified close to PsCam037030_22140_221, one of which being also precisely annotated as a *CBF14* by Tayeh et al. [36]. In the LD block VI.2, 2 FD-associated markers, namely PsCam023246_13111_1125 and PsCam007060_5248_2156, were also found to be close to AP2/ERF (Ethylene-responsive transcription factor) genes, namely Psat1g097280 and Psat1g097280, for which a precise study of the sequence is needed to verify if they belong to the CBF sub-family. Finally, 4 other potential candidate genes were found at the vicinity of the LD block VI.2, namely Psat1g103920, Psat1g106640, Psat1g115640 and Psat1g103680. These genes were annotated as *DREB* (Dehydration-responsive element-binding protein) genes, another term used to refer to *CBF* genes, in the homologous genes of *M. truncatula* (Additional file 11). These 4 genes lie in an interval of 35 Mb situated at a distance of 448 kb from the three other *CBF* genes mentioned above. One of them, Psat1g103920, contained the marker PsCam050192_32788_145 which was excluded from the GWA analysis because it showed a minor allele frequency lower than 0.05.

Three FD-associated markers, namely PsCam035617_20792_637 (LD-block III.1), PsCam048068_30823_2326 (LD-block V.1) and PsCam011774_8038_200 (LD-block VI.8) were found to correspond to 3 genes, namely Psat5g299600, Psat3g087400 and Psat1g119400 encoding brassinosteroid receptors (Additional file 11).

The FD-related marker PsCam037922_22979_691, which is only 0.1 cM apart from PsCam035617_20792_637 just mentioned above, was found to lie in Psat5g299720, a gene encoding for the gibberellin 3beta-hydroxylase enzyme and shown to correspond to the dwarfism gene *Le* in pea (Additional file 11).

Finally, three FD-associated markers belonging to the FD-related locus on LG VII, namely PsCam001108_940_48, PsCam037927_22984_97 and PsCam004928_3732_3087 were located within genes related to the synthesis of soluble sugars (Psat7g180280: endo-beta-1,3-glucanase, Psat7g193120: endo-1,4-beta-glucanase and Psat7g214880: beta-glucosidase G2, respectively) (Additional file 11).

## Discussion

**GWAS brings new insights into the determinism of frost tolerance in pea**

In the present study, 75 markers associated with frost tolerance, 62 markers significantly detected by GWAS

and 13 markers in high LD ($r^2 > 0.8$) with one or the other of the 62 markers, were located among all the linkage groups of pea, except LG IV (Table 3, Additional file 9).

Comparing the map positions with those of frost tolerance QTLs previously described, 3 regions corresponding to 62 markers among the 75 markers, were found to colocalize with 3 main QTLs previously detected by linkage mapping in two bi-parental populations for the same trait, namely WFD 3.2, WFD 5.1 and WFD 6.1 [21–23] (Fig. 4). These 3 QTLs were repeatedly detected in 5, 11 and 10 field conditions for WFD 3.2, WFD 5.1 and WFD 6.1, respectively. Moreover, the position corresponding to WFD 6.1 seems to also match with an EST marker recently found to be associated with frost tolerance in pea [38]: indeed, Liu et al. identified 7 marker-trait associations within a collection of 672 accessions, among which the marker EST1109 was located on LG VI within a functional gene that has a high homology with a gene encoding an alpha-mannosidase in *M. truncatula*. We reviewed the 1646 transcript-derived SNP markers mapped on LG VI by Tayeh et al. [36] and found that a marker corresponding to an alpha-mannosidase-like protein was located at 43.2 cM on the consensus map; consequently, as the FD-associated markers detected on LG VI in the present study are comprised between 41.4 and 52.8 cM, it is likely that the LG VI positions of both association studies coincide. Similarly, this study validated the QTL previously detected on LG III (WFD 3.2 in Pop2 [21], III.2 in Pop9 [23]) with a higher resolution than previous linkage mapping studies by the identifications of three main significant FD-associated markers, and whose positive alleles decreased the frost damage of accessions by an average of 0.33 (Table 3). All favorable alleles were carried by the sensitive genotypes Térèse and Caméor unlike the tolerant genotypes Champagne and China which had all undesirable alleles underlying this locus (Table 3). Altogether, the consistency of the 3 positions detected in both bi-parental mapping and association genetics reinforces their interest for breeding. The results presented here constitute an additional step towards the identification of underlying genes potentially involved in the control of frost tolerance, thanks to refined intervals provided by GWAS.

Unlike the correspondences with bi-parental mapping positions presented above, the present GWA study did not highlight any colocalization with a major QTL, namely WFD3.1, which was however found to be responsible for up to 52 and 19% of the winter frost damage variation within the RILs populations derived from Champagne x Térèse (Pop2) and China x Caméor (Pop9), respectively [21, 23]. The flowering gene *Hr* (*High response to photoperiod*), an orthologue of the Arabidopsis *Early flowering 3*, *Elf3* [40], was shown to be
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Comparative genetic map of genome-wide association study (GWAS) loci identified in the present study and quantitative trait loci (QTL) previously detected for frost tolerance in pea. Only linkage groups (LGs) with significant frost damages (FD)-associated markers detected by GWAS are presented. On each LG, SNP markers are shown on the right and genetic positions between markers are indicated in cM on the left. Frost tolerance loci detected in the present study are shown in red: FD-associated markers are in a red and underlined font; markers in linkage disequilibrium (LD; r² > 0.8) with associated marker(s) are in a red and non-underlined font; the LD blocks identified by GWAS are drawn as red bars on the right of each LG. QTLs represented by blue and green bars were detected in the Champagne x Térèse [21, 22] and China x Caméor [23] populations respectively. For presentation purposes, only markers at the vicinity of significant loci and a few markers distributed along LGs are shown</figcaption></figure>

a relevant candidate for this QTL as it allows plants to be maintained in a vegetative state under short days and thus to escape the main winter freezing periods. It is likely that, in the case of the WFD3.1 position corresponding to the *Hr* gene, a strong correlation may have emerged between the population structure, possibly biased by the allelic variation at the *Hr* locus and the frost damage trait. This hypothesis relies on the observation that *Hr* may have been the target of natural selection for frost tolerance. Weller et al. [40] speculated that the *hr* mutation may have arisen within an ancestral pea lineage originating from the Near East domestication center and carrying the *Hr* allele. The *hr* mutation possibly permitted summer cropping in areas characterized by colder winters and is therefore highly represented in many domesticated lines of *Pisum sativum* at the origin of the current spring peas. To explore the hypothesis of undetection of a true association due to confounding with the population structure, we have verified the distribution of the *Hr* alleles, represented by their *Elf3* genotype, within the DAPC clusters of the association mapping population (Additional file 12). The *Hr* accessions, homozygous for the dominant allele *Hr*, are the main components of the clusters 7 (96%) and 1 (84%) and they represent 56.9% of the cluster 4. They are thus over represented in the three clusters gathering most of the winter sowing-type accessions, which may have contributed to a correlation between the frost tolerance trait and the population structure. The *hr* accessions, homozygous for the recessive allele *hr*, are the main components of the clusters 3 (100%), 2 (100%), 5 (98.7%) and 6 (95.7%), which are mostly spring sowing-type accessions. Consequently, we can suggest that the correction for the population structure (Q matrix) might have resulted in a structuring marker that probably cannot be detected by further association analysis using this sample of accessions. This kind of result has already been reported by Visioni et al. [27] who found that 2 of the 3 most significant SNPs of their study, tightly linked to major known
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Comparative genetic map of genome-wide association study (GWAS) loci identified in the present study and quantitative trait loci (QTL) previously detected for frost tolerance in pea. Only linkage groups (LGs) with significant frost damages (FD)-associated markers detected by GWAS are presented. On each LG, SNP markers are shown on the right and genetic positions between markers are indicated in cM on the left. Frost tolerance loci detected in the present study are shown in red: FD-associated markers are in a red and underlined font; markers in linkage disequilibrium (LD; r² > 0.8) with associated marker(s) are in a red and non-underlined font; the LD blocks identified by GWAS are drawn as red bars on the right of each LG. QTLs represented by blue and green bars were detected in the Champagne x Térèse [21, 22] and China x Caméor [23] populations respectively. For presentation purposes, only markers at the vicinity of significant loci and a few markers distributed along LGs are shown</figcaption></figure>

a relevant candidate for this QTL as it allows plants to be maintained in a vegetative state under short days and thus to escape the main winter freezing periods. It is likely that, in the case of the WFD3.1 position corresponding to the *Hr* gene, a strong correlation may have emerged between the population structure, possibly biased by the allelic variation at the *Hr* locus and the frost damage trait. This hypothesis relies on the observation that *Hr* may have been the target of natural selection for frost tolerance. Weller et al. [40] speculated that the *hr* mutation may have arisen within an ancestral pea lineage originating from the Near East domestication center and carrying the *Hr* allele. The *hr* mutation possibly permitted summer cropping in areas characterized by colder winters and is therefore highly represented in many domesticated lines of *Pisum sativum* at the origin of the current spring peas. To explore the hypothesis of undetection of a true association due to confounding with the population structure, we have verified the distribution of the *Hr* alleles, represented by their *Elf3* genotype, within the DAPC clusters of the association mapping population (Additional file 12). The *Hr* accessions, homozygous for the dominant allele *Hr*, are the main components of the clusters 7 (96%) and 1 (84%) and they represent 56.9% of the cluster 4. They are thus over represented in the three clusters gathering most of the winter sowing-type accessions, which may have contributed to a correlation between the frost tolerance trait and the population structure. The *hr* accessions, homozygous for the recessive allele *hr*, are the main components of the clusters 3 (100%), 2 (100%), 5 (98.7%) and 6 (95.7%), which are mostly spring sowing-type accessions. Consequently, we can suggest that the correction for the population structure (Q matrix) might have resulted in a structuring marker that probably cannot be detected by further association analysis using this sample of accessions. This kind of result has already been reported by Visioni et al. [27] who found that 2 of the 3 most significant SNPs of their study, tightly linked to major known
genetic determinants of cold tolerance in barley, were
undetected by GWAS if a correction by the structure was
used. It was particularly the case for a SNP linked to *Vrn-H1*, a developmental locus governing barley vernalization
requirement, which is for long a candidate for the frost
tolerance locus *Fr-H1* but whose effect was suspected to be
confounded with the population structure. To overcome
this point, NAM-like linkage populations with bi-parental
crosses in a reference design could be an interesting plant
material for association mapping, in order to minimize
population structure which may be necessary for dissecting
the most structured traits.

Comparatively with consistent QTLs previously de-
tected in bi-parental populations, the present GWA study
also pinpointed three loci which have either not been de-
tected yet (one region on LG II) or formerly detected in
only one environment (two regions located on LG I and
LG VII, respectively). The significant marker on LG II is
supported only by one experiment in controlled condi-
tions and must therefore be used with caution in breeding
programs even if two distinct favorable and unfavorable
haplotypes were identified (Additional file 10). The two
markers significantly associated with frost tolerance on
LG I are located at 47.5 cM on the consensus map, which
lies in the projected confidence interval of a previous
QTL, namely WFDcle.a, identified in one field condition
with the *Hr* subpopulation extracted from Pop2 [22]. This
colocalization slightly reinforces the consistency of this
LG I position. In the same way, the 8 significant markers
identified on LG VII with this study were located between
72.9 and 89.3 cM on the consensus map, which overlap
with one former QTL, namely FD.c, which was detected
once in a controlled chamber experiment [22], with the
same *Hr* subpopulation. Thus, the colocalizing region on
LG VII relies now on three independent experiments. In a
panel of 672 accessions, Liu et al. [38] identified one
marker on LG I and two markers on LG VII that were sig-
nificantly associated with frost tolerance. It would be in-
teresting to check the localization of these markers on the
consensus map used here, as their position on the genetic
map used by the authors [41] is not reported. Additionally,
to these coincident positions for the frost tolerance trait
on LG I and LG VII within different experimental condi-
tions, the LG I and the LG VII regions detected in this
study seems to also overlap with already detected QTLs
for resistance to *Aphanomyces euteiches*. Indeed, both
markers of the FD-related LD block I.1 lie in the confi-
dence interval of the QTL *Ae-Ps1.1* identified by Hamon
et al. [42], when the latter is projected on the consensus
map used in the present work. Besides, the FD-related
locus on LG VII includes the LD block VII.16 associated
to the resistance to *Aphanomyces euteiches* [35], with
which it shares the FD-associated marker PsCam038378_
23415_721. The FD-related locus on LG VII deserves a

particular attention for a further use in breeding for frost
tolerance because accessions carrying the favorable haplo-
type underlying this locus (haplotype VII.4) shows very
low scores of frost damages, ranging from 0 to 0.96,
both in the field and controlled conditions experiments
(Additional file 10). The relationships between haplo-
types at this locus and the values for frost tolerance
and *Aphanomyces euteiches* resistance will however
have to be more precisely explored, to check if both
traits can be bred favorably at this locus.

GWAS detected frost tolerance-associated markers which
are included in relevant putative candidate genes

The projection of the 75 FD-related markers on the pea
genome assembly [39] allowed to identify 590 annotated
genes with known putative protein functions, located in
an interval of ±1 Mb on both sides of FD-related markers
(Additional file 11). Among the diverse protein functions
predicted, some are already related to the acquisition of
frost tolerance in the literature.

Comparison of map positions has shown that the FD-
related locus detected on LG VI in this study colocalizes
with the previous QTL WFD 6.1, which is itself ortholo-
gous to a major QTL for frost tolerance in *M. trunca-
tula* (Mt-FTQTL6) [43]. Tayeh et al. moreover showed
that Mt-FTQTL6 covers a region containing a cluster of
twelve *CBF* genes tandemly duplicated [19]. In the
present study, 9 *AP2* domain genes were found to cor-
respond, or to be at the vicinity of, FD-related markers
of LG VI. Among these genes, 7 are annotated as *CBF*
or *DREB* genes. Given these results and the previous
findings of Tayeh et al. concerning Mt-FTQTL6 [19],
*CBF* genes located in the LD block VI.2, or at its vicinity,
are also relevant candidates determining frost tolerance
at this locus in pea. The potential role of *CBF* genes, and
particularly *CBF14*, has already been highlighted in
cereals. In wheat, Zhu et al. [44] showed that the natural
variation for frost tolerance is mainly associated with a
*frost resistance 2* (*FR2*) locus including tandemly repli-
cated *CBF* genes that regulates the expression of cold-
regulated genes. Additionally, these authors proved that
an increased copy number of *CBF14* was frequently as-
sociated with the tolerant haplotype of the locus *FR-A2*
and with higher *CBF14* transcript levels in response to
cold. Novák et al. [45] showed that *CBF14* genes con-
tribute to enhance frost tolerance during cold acclima-
tion in cereals.

Three candidate genes corresponding to FD-associated
markers detected on LG III, LG V and LG VI, and anno-
tated as brassinosteroid receptors, also appear in the
literature to be implied in the crosstalk between plant
hormone signaling in the cold stress response and the
CBF regulon. In Arabidopsis, Eremina et al. [46] pro-
vided evidence that brassinosteroids contribute to the
genetic determinants of cold tolerance in barley, were
undetected by GWAS if a correction by the structure was
used. It was particularly the case for a SNP linked to *Vrn-H1*, a developmental locus governing barley vernalization
requirement, which is for long a candidate for the frost
tolerance locus *Fr-H1* but whose effect was suspected to be
confounded with the population structure. To overcome
this point, NAM-like linkage populations with bi-parental
crosses in a reference design could be an interesting plant
material for association mapping, in order to minimize
population structure which may be necessary for dissecting
the most structured traits.

Comparatively with consistent QTLs previously de-
tected in bi-parental populations, the present GWA study
also pinpointed three loci which have either not been de-
tected yet (one region on LG II) or formerly detected in
only one environment (two regions located on LG I and
LG VII, respectively). The significant marker on LG II is
supported only by one experiment in controlled condi-
tions and must therefore be used with caution in breeding
programs even if two distinct favorable and unfavorable
haplotypes were identified (Additional file 10). The two
markers significantly associated with frost tolerance on
LG I are located at 47.5 cM on the consensus map, which
lies in the projected confidence interval of a previous
QTL, namely WFDcle.a, identified in one field condition
with the *Hr* subpopulation extracted from Pop2 [22]. This
colocalization slightly reinforces the consistency of this
LG I position. In the same way, the 8 significant markers
identified on LG VII with this study were located between
72.9 and 89.3 cM on the consensus map, which overlap
with one former QTL, namely FD.c, which was detected
once in a controlled chamber experiment [22], with the
same *Hr* subpopulation. Thus, the colocalizing region on
LG VII relies now on three independent experiments. In a
panel of 672 accessions, Liu et al. [38] identified one
marker on LG I and two markers on LG VII that were sig-
nificantly associated with frost tolerance. It would be in-
teresting to check the localization of these markers on the
consensus map used here, as their position on the genetic
map used by the authors [41] is not reported. Additionally,
to these coincident positions for the frost tolerance trait
on LG I and LG VII within different experimental condi-
tions, the LG I and the LG VII regions detected in this
study seems to also overlap with already detected QTLs
for resistance to *Aphanomyces euteiches*. Indeed, both
markers of the FD-related LD block I.1 lie in the confi-
dence interval of the QTL *Ae-Ps1.1* identified by Hamon
et al. [42], when the latter is projected on the consensus
map used in the present work. Besides, the FD-related
locus on LG VII includes the LD block VII.16 associated
to the resistance to *Aphanomyces euteiches* [35], with
which it shares the FD-associated marker PsCam038378_
23415_721. The FD-related locus on LG VII deserves a

particular attention for a further use in breeding for frost
tolerance because accessions carrying the favorable haplo-
type underlying this locus (haplotype VII.4) shows very
low scores of frost damages, ranging from 0 to 0.96,
both in the field and controlled conditions experiments
(Additional file 10). The relationships between haplo-
types at this locus and the values for frost tolerance
and *Aphanomyces euteiches* resistance will however
have to be more precisely explored, to check if both
traits can be bred favorably at this locus.

GWAS detected frost tolerance-associated markers which
are included in relevant putative candidate genes

The projection of the 75 FD-related markers on the pea
genome assembly [39] allowed to identify 590 annotated
genes with known putative protein functions, located in
an interval of ±1 Mb on both sides of FD-related markers
(Additional file 11). Among the diverse protein functions
predicted, some are already related to the acquisition of
frost tolerance in the literature.

Comparison of map positions has shown that the FD-
related locus detected on LG VI in this study colocalizes
with the previous QTL WFD 6.1, which is itself ortholo-
gous to a major QTL for frost tolerance in *M. trunca-
tula* (Mt-FTQTL6) [43]. Tayeh et al. moreover showed
that Mt-FTQTL6 covers a region containing a cluster of
twelve *CBF* genes tandemly duplicated [19]. In the
present study, 9 *AP2* domain genes were found to cor-
respond, or to be at the vicinity of, FD-related markers
of LG VI. Among these genes, 7 are annotated as *CBF*
or *DREB* genes. Given these results and the previous
findings of Tayeh et al. concerning Mt-FTQTL6 [19],
*CBF* genes located in the LD block VI.2, or at its vicinity,
are also relevant candidates determining frost tolerance
at this locus in pea. The potential role of *CBF* genes, and
particularly *CBF14*, has already been highlighted in
cereals. In wheat, Zhu et al. [44] showed that the natural
variation for frost tolerance is mainly associated with a
*frost resistance 2 (FR2)* locus including tandemly repli-
cated *CBF* genes that regulates the expression of cold-
regulated genes. Additionally, these authors proved that
an increased copy number of *CBF14* was frequently as-
sociated with the tolerant haplotype of the locus *FR-A2*
and with higher *CBF14* transcript levels in response to
cold. Novák et al. [45] showed that *CBF14* genes con-
tribute to enhance frost tolerance during cold acclima-
tion in cereals.

Three candidate genes corresponding to FD-associated
markers detected on LG III, LG V and LG VI, and anno-
tated as brassinosteroid receptors, also appear in the
literature to be implied in the crosstalk between plant
hormone signaling in the cold stress response and the
CBF regulon. In Arabidopsis, Eremina et al. [46] pro-
vided evidence that brassinosteroids contribute to the
control of freezing tolerance. Indeed, these authors showed that brassinosteroid-deficient mutants of Arabidopsis were hypersensitive to freezing stress, whereas an activation of the brassinosteroid signaling pathway increased freezing tolerance both before and after cold acclimation. Furthermore, two brassinosteroids-responsive transcription factors have also been characterized as direct regulators of CBF expression through their binding to the promoters of these genes [46, 47]. In cultivated plants, the role of brassinosteroids is so far documented for the response to the chilling stress, as reviewed by Anwar et al. [48].

Within the FD-related locus of LG III, was identified a candidate gene encoding for the gibberellin 3beta-hydroxylase enzyme which produces bioactive gibberellin, also known as Le in pea (Additional file 11). Recessive *le* mutants at this locus are impaired in the production of gibberellin and produce a dwarf phenotype [49, 50]. In Arabidopsis, Achard et al. [51] undertook a molecular and genetic approach to evaluate the interaction between the *CBF1*-dependent cold acclimation pathway and the gibberellin pathway. They proposed a model in which the induction of *CBF1* expression by low temperature affects the gibberellin metabolism via upregulation of gibberellin-2-oxydase gene transcripts. The following reduction in bioactive gibberellin causes a higher accumulation of DELLAs, a family of nuclear growth-repressing proteins which in turns restrains plant growth. The *Le* gene has already been proposed as a candidate for the WFD3.2 and III.3 QTLs identified in the Pop2 and Pop9 populations, respectively. We considered the haplotypes of the parents of both populations at the corresponding FD-associated locus in this study and found that the favorable haplotype (III.1, including the dwarf allele at the *Le* gene) was borne by Térèse and Caméor, while the unfavorable haplotype (III.2, including the wild-type allele at the *Le* gene) was carried by Champagne and China. This observation is consistent with the favorable and unfavorable alleles determined by QTL mapping in bi-parental populations. As the three genes constituting the FD-related locus of LG III lie within a 1 cM interval, neither a synergistic effect nor linkage can be excluded.

Finally, three candidates underlying the FD-related locus on LG VII (chromosome 7) corresponded to genes related to the synthesis of soluble sugars. As accumulation of soluble sugars during cold acclimation is well documented in many plants, we can suggest that the LG VII-locus may have a role in frost tolerance by accumulating sugars in plant tissues during cold acclimation. Several roles for sugars in protecting cells from freezing injury have been proposed, including functioning as cryoprotectants for specific enzymes, as molecules promoting membrane stability and as osmolytes to prevent excessive dehydration during freezing, as reviewed by Xin and Browse [52].

## GWAS provides new markers and new genitors to breed for frost tolerance in pea

In the present study, 6 loci related to frost tolerance in pea were identified. At these 6 FD-related loci, 7 favorable haplotypes, carrying the highest number of favorable alleles at the FD-associated markers detected by GWAS, were significantly associated with the lowest scores of frost damages, i.e. the highest levels of frost tolerance. We identified 12 accessions showing lower scores of mean frost damages ranging from 0.13 to 1.04 and cumulating 6 (Glacier), 5 (Melrose, Blixt 7, Winterberger, Holly 9, Black seeded, Holly 17, Blixt 109, Fe and P1259) or 4 (Cote d'or and Picar) favorable haplotypes at the 6 FD-related loci (Additional file 13). All these accessions belong to the cluster 7 which was shown to be totally isolated from the other clusters by the DAPC analysis. They are fodder peas among which frost tolerant accessions have already been identified for breeding. The same 12 accessions also carry the *Hr* allele which was formerly shown to be favorable to frost tolerance [21]. One common point of these accessions except the accession named 'Glacier, is that they carry the unfavorable haplotype at the FD-related locus on LG III comprising the *Le* gene. But rather than indicating a minor effect of the favorable haplotype at this locus, genetically linked to the dwarf *le* allele, this feature is to relate to the observation that field-autumn-sown *Hr* lines remain dwarf until a longer spring daylength has also triggered off the switch from the prostrate to the erected growth habit. This suggests an epistatic effect of *Hr* upon the expression of the dwarfism [21]. Comparatively to the above-described material, we also identified 16 accessions carrying all the unfavorable alleles at the FD-related loci located on LG V, VI and VII of the present study, additionally to the unfavorable allele *hr*. These accessions presented only 3 (Petit provençal, eM and Cador), 2 (Pi196033, Aldot, Chine-d368, 667, Merveille de Kelvedon, Miravil, Wav f502, Mingomark) or 1 (Ceia, Alaska, Finlande, Ersling and Automobile) favorable haplotype(s) identified on LG I, II and/or III, and which are mainly garden spring cultivars or breeding accessions presenting higher scores of mean frost damages ranging from 2.38 and 4.56. This allows us to state that the three loci on LG V, VI and VII play a bigger part in the frost tolerance on pea than the other FD-related loci located on LG I, II and III. Our results can help choosing tolerant progenitors and following favorable haplotypes through marker-assisted breeding. Furthermore, the FD-associated locus on LGV was found to overlap with the confidence interval of the frost tolerance QTL WFD 5.1 earlier detected within the Pop2 population [21] (Fig. 4). Comparatively to the linkage mapping method, with which a confidence interval of 16.6 cM was obtained, the GWA study enabled to refine the confidence interval to
control of freezing tolerance. Indeed, these authors showed that brassinosteroid-deficient mutants of Arabidopsis were hypersensitive to freezing stress, whereas an activation of the brassinosteroid signaling pathway increased freezing tolerance both before and after cold acclimation. Furthermore, two brassinosteroids-responsive transcription factors have also been characterized as direct regulators of CBF expression through their binding to the promoters of these genes [46, 47]. In cultivated plants, the role of brassinosteroids is so far documented for the response to the chilling stress, as reviewed by Anwar et al. [48].

Within the FD-related locus of LG III, was identified a candidate gene encoding for the gibberellin 3beta-hydroxylase enzyme which produces bioactive gibberellin, also known as Le in pea (Additional file 11). Recessive *le* mutants at this locus are impaired in the production of gibberellin and produce a dwarf phenotype [49, 50]. In Arabidopsis, Achard et al. [51] undertook a molecular and genetic approach to evaluate the interaction between the *CBF1*-dependent cold acclimation pathway and the gibberellin pathway. They proposed a model in which the induction of *CBF1* expression by low temperature affects the gibberellin metabolism via upregulation of gibberellin-2-oxydase gene transcripts. The following reduction in bioactive gibberellin causes a higher accumulation of DELLAs, a family of nuclear growth-repressing proteins which in turns restrains plant growth. The *Le* gene has already been proposed as a candidate for the WFD3.2 and III.3 QTLs identified in the Pop2 and Pop9 populations, respectively. We considered the haplotypes of the parents of both populations at the corresponding FD-associated locus in this study and found that the favorable haplotype (III.1, including the dwarf allele at the *Le* gene) was borne by Térèse and Caméor, while the unfavorable haplotype (III.2, including the wild-type allele at the *Le* gene) was carried by Champagne and China. This observation is consistent with the favorable and unfavorable alleles determined by QTL mapping in bi-parental populations. As the three genes constituting the FD-related locus of LG III lie within a 1 cM interval, neither a synergistic effect nor linkage can be excluded.

Finally, three candidates underlying the FD-related locus on LG VII (chromosome 7) corresponded to genes related to the synthesis of soluble sugars. As accumulation of soluble sugars during cold acclimation is well documented in many plants, we can suggest that the LG VII-locus may have a role in frost tolerance by accumulating sugars in plant tissues during cold acclimation. Several roles for sugars in protecting cells from freezing injury have been proposed, including functioning as cryoprotectants for specific enzymes, as molecules promoting membrane stability and as osmolytes to prevent excessive dehydration during freezing, as reviewed by Xin and Browse [52].

GWAS provides new markers and new genitors to breed for frost tolerance in pea

In the present study, 6 loci related to frost tolerance in pea were identified. At these 6 FD-related loci, 7 favorable haplotypes, carrying the highest number of favorable alleles at the FD-associated markers detected by GWAS, were significantly associated with the lowest scores of frost damages, i.e. the highest levels of frost tolerance. We identified 12 accessions showing lower scores of mean frost damages ranging from 0.13 to 1.04 and cumulating 6 (Glacier), 5 (Melrose, Blixt 7, Winterberger, Holly 9, Black seeded, Holly 17, Blixt 109, Fe and P1259) or 4 (Cote d'or and Picar) favorable haplotypes at the 6 FD-related loci (Additional file 13). All these accessions belong to the cluster 7 which was shown to be totally isolated from the other clusters by the DAPC analysis. They are fodder peas among which frost tolerant accessions have already been identified for breeding. The same 12 accessions also carry the *Hr* allele which was formerly shown to be favorable to frost tolerance [21]. One common point of these accessions except the accession named 'Glacier, is that they carry the unfavorable haplotype at the FD-related locus on LG III comprising the *Le* gene. But rather than indicating a minor effect of the favorable haplotype at this locus, genetically linked to the dwarf *le* allele, this feature is to relate to the observation that field-autumn-sown *Hr* lines remain dwarf until a longer spring daylength has also triggered off the switch from the prostrate to the erected growth habit. This suggests an epistatic effect of *Hr* upon the expression of the dwarfism [21]. Comparatively to the above-described material, we also identified 16 accessions carrying all the unfavorable alleles at the FD-related loci located on LG V, VI and VII of the present study, additionally to the unfavorable allele *hr*. These accessions presented only 3 (Petit provençal, eM and Cador), 2 (Pi196033, Aldot, Chine-d368, 667, Merveille de Kelvedon, Miravil, Wav f502, Mingomark) or 1 (Ceia, Alaska, Finlande, Ersling and Automobile) favorable haplotype(s) identified on LG I, II and/or III, and which are mainly garden spring cultivars or breeding accessions presenting higher scores of mean frost damages ranging from 2.38 and 4.56. This allows us to state that the three loci on LG V, VI and VII play a bigger part in the frost tolerance on pea than the other FD-related loci located on LG I, II and III. Our results can help choosing tolerant progenitors and following favorable haplotypes through marker-assisted breeding. Furthermore, the FD-associated locus on LGV was found to overlap with the confidence interval of the frost tolerance QTL WFD 5.1 earlier detected within the Pop2 population [21] (Fig. 4). Comparatively to the linkage mapping method, with which a confidence interval of 16.6 cM was obtained, the GWA study enabled to refine the confidence interval to
7.4 cM (53.6 to 61 cM on the consensus map). This refined LG V locus presents a particular interest in breeding as it may provide markers to break the genetic relationship between the frost tolerance position and the neighbouring locus governing the seed trypsin activity. Trypsin inhibitors are known to be unfavorable for animal feed because they decrease the digestibility of protein [53]. The locus responsible for the seed trypsin activity (Tri) has been mapped on LG V [54] within the confidence interval of WFD 5.1 [21]. The favorable alleles at both loci are generally in repulsion. On the consensus map, this locus is represented by three markers, annotated as *trypsin inhibitor* genes [36] and located between 67.0 and 67.3 cM, 6 cM apart from the frost tolerance locus detected by GWAS. Thus, it seems possible to select favorable alleles for frost tolerance corresponding to the FD-associated markers detected by GWAS on LG V together with recessive alleles of the markers encoding for the trypsin inhibitor genes.

## Conclusion

In the present study, GWAS enabled to confirm QTLs significantly associated with frost tolerance such as WFD 3.2, WFD 5.1 and WFD 6.1. It also allowed to identify one region on LG II, which has not been detected yet and provided significant associations for two regions on LGI and LG VII that were formerly detected in only one environment. The results showed that GWAS is an effective strategy to identify markers precisely defining frost tolerance loci, which can be useful to breed for antagonistic traits as it is for the frost tolerance and Tri loci on LG V which are in linkage disequilibrium and in a repulsion phase. Our results also highlight that GWAS enables to find new sources of frost tolerance within collections of pea genetic resources. Finally, the present GWA study also brought to light the presence of CBF transcriptions factors as potential genetic determinants of the frost tolerance locus on LG VI, with one CBF-annotated marker being in high LD with significant FD-associated markers of the locus and six additional *CBF/DREB*-annotated genes mapped at the vicinity. As 12 tandemly duplicated *CBF* genes were already found to be relevant candidates underlying the orthologous frost tolerance QTL on *Medicago truncatula* chromosome 6, the hypothesis of a similar genomic organization in pea deserves to be tested.

## Methods

### Plant material

The association mapping collection, also named the collection, consists of 365 accessions (Additional file 6) from the pea reference collection described in Burstin et al. [55]. *Pisum* accessions from the collection represents a large genetic diversity ranging from wild peas

(*Pisum fulvum*, *P. humile*, *P. elatius*, *P. speciosum*, *P. transcaucasicum* and *P. abyssinicum*) and landraces, to breeding lines and cultivars. This collection also represents a variability of genotypes based on the type of sowing (winter vs spring peas) and the type of end-use (fodder, field, mangetout, preserve and garden peas). Reference accessions which are the parents of bi-parental populations formerly used in QTL studies for frost tolerance [21–23], i.e. Champagne, China, Térèse and Caméor, are included in the collection.

All genotypes were purified for one generation by single seed descent (SSD) in insect-proof glasshouse. After this SSD generation, seeds were increased for one generation in insect-proof glasshouse. The seeds produced were sown in a nursery: tissue samples were harvested in bulk for DNA production from 10 sister plants and harvested seeds were used for phenotyping. When necessary, DNA was extracted again from the offspring of these plants. There is therefore zero or one generation between phenotyping and genotyping.

### Phenotyping

Frost tolerance of the collection was evaluated under field and controlled conditions. The field experiment was carried out at the INRAE (National Research Institute for Agriculture, Food and Environment) experimental station of Clermont-Ferrand Theix, France (45.72°N latitude and 3.02°E longitude at an altitude of 890 m) during the growing season of 2007–2008. Sowing date was 09 October 2007 and the date of emergence was 26 October 2007. Plots were sown in a randomized complete block design with three replicates. Weeds and diseases were controlled chemically.

The record of temperatures indicated that cold acclimation and freezing periods occurred during the experiment (Additional file 14). The collection was assessed for frost tolerance by visual estimation of winter frost damages after the main winter freezing periods have passed. As described in previous studies [21–23], a score was assigned to a plot as a whole, based on the extent of necrotic areas of the aerial parts of the plants. The scale ranged from 0 to 5 where 0 represented no damage and 5 a dead plant. Frost damages observations were realized at four dates in 2008: January, 4th and 15th, March, 28th and April, 10th.

A frost experiment was also conducted in a controlled environment chamber using the standardized test described previously by Dumont et al. [22], which mimics the successive periods of cold acclimation and frost generally encountered in the field by autumn-sown peas. Pea accessions were placed according to a randomized complete block design with three replicates. To provide three biological replicates, the experiment was carried out three times successively in the same controlled
7.4 cM (53.6 to 61 cM on the consensus map). This refined LG V locus presents a particular interest in breeding as it may provide markers to break the genetic relationship between the frost tolerance position and the neighbouring locus governing the seed trypsin activity. Trypsin inhibitors are known to be unfavorable for animal feed because they decrease the digestibility of protein [53]. The locus responsible for the seed trypsin activity (Tri) has been mapped on LG V [54] within the confidence interval of WFD 5.1 [21]. The favorable alleles at both loci are generally in repulsion. On the consensus map, this locus is represented by three markers, annotated as *trypsin inhibitor* genes [36] and located between 67.0 and 67.3 cM, 6 cM apart from the frost tolerance locus detected by GWAS. Thus, it seems possible to select favorable alleles for frost tolerance corresponding to the FD-associated markers detected by GWAS on LG V together with recessive alleles of the markers encoding for the trypsin inhibitor genes.

## Conclusion

In the present study, GWAS enabled to confirm QTLs significantly associated with frost tolerance such as WFD 3.2, WFD 5.1 and WFD 6.1. It also allowed to identify one region on LG II, which has not been detected yet and provided significant associations for two regions on LGI and LG VII that were formerly detected in only one environment. The results showed that GWAS is an effective strategy to identify markers precisely defining frost tolerance loci, which can be useful to breed for antagonistic traits as it is for the frost tolerance and Tri loci on LG V which are in linkage disequilibrium and in a repulsion phase. Our results also highlight that GWAS enables to find new sources of frost tolerance within collections of pea genetic resources. Finally, the present GWA study also brought to light the presence of CBF transcriptions factors as potential genetic determinants of the frost tolerance locus on LG VI, with one CBF-annotated marker being in high LD with significant FD-associated markers of the locus and six additional *CBF/DREB*-annotated genes mapped at the vicinity. As 12 tandemly duplicated *CBF* genes were already found to be relevant candidates underlying the orthologous frost tolerance QTL on *Medicago truncatula* chromosome 6, the hypothesis of a similar genomic organization in pea deserves to be tested.

## Methods

### Plant material

The association mapping collection, also named the collection, consists of 365 accessions (Additional file 6) from the pea reference collection described in Burstin et al. [55]. *Pisum* accessions from the collection represents a large genetic diversity ranging from wild peas

(*Pisum fulvum*, *P. humile*, *P. elatius*, *P. speciosum*, *P. transcaucasicum* and *P. abyssinicum*) and landraces, to breeding lines and cultivars. This collection also represents a variability of genotypes based on the type of sowing (winter vs spring peas) and the type of end-use (fodder, field, mangetout, preserve and garden peas). Reference accessions which are the parents of bi-parental populations formerly used in QTL studies for frost tolerance [21–23], i.e. Champagne, China, Térèse and Caméor, are included in the collection.

All genotypes were purified for one generation by single seed descent (SSD) in insect-proof glasshouse. After this SSD generation, seeds were increased for one generation in insect-proof glasshouse. The seeds produced were sown in a nursery: tissue samples were harvested in bulk for DNA production from 10 sister plants and harvested seeds were used for phenotyping. When necessary, DNA was extracted again from the offspring of these plants. There is therefore zero or one generation between phenotyping and genotyping.

### Phenotyping

Frost tolerance of the collection was evaluated under field and controlled conditions. The field experiment was carried out at the INRAE (National Research Institute for Agriculture, Food and Environment) experimental station of Clermont-Ferrand Theix, France (45.72°N latitude and 3.02°E longitude at an altitude of 890 m) during the growing season of 2007–2008. Sowing date was 09 October 2007 and the date of emergence was 26 October 2007. Plots were sown in a randomized complete block design with three replicates. Weeds and diseases were controlled chemically.

The record of temperatures indicated that cold acclimation and freezing periods occurred during the experiment (Additional file 14). The collection was assessed for frost tolerance by visual estimation of winter frost damages after the main winter freezing periods have passed. As described in previous studies [21–23], a score was assigned to a plot as a whole, based on the extent of necrotic areas of the aerial parts of the plants. The scale ranged from 0 to 5 where 0 represented no damage and 5 a dead plant. Frost damages observations were realized at four dates in 2008: January, 4th and 15th, March, 28th and April, 10th.

A frost experiment was also conducted in a controlled environment chamber using the standardized test described previously by Dumont et al. [22], which mimics the successive periods of cold acclimation and frost generally encountered in the field by autumn-sown peas. Pea accessions were placed according to a randomized complete block design with three replicates. To provide three biological replicates, the experiment was carried out three times successively in the same controlled
environment chamber. The temperature, light level and humidity were recorded and were similar during the three experiments. Briefly, the plants at the stage of 2nd - 3rd leaf were first treated with a regime of 11 days of cold acclimation at 10 °C/2 °C (day/night) with a 10 h photoperiod. The frost treatment was then carried out at 6 °C/- 8 °C with 8 h of daylight during 4 days. After frost, a recovery period was applied with a temperature regime of 16 °C/5 °C and 10 h of daylight during 8 days. Frost tolerance was evaluated by scoring frost damages at the end of the recovery period with the same scoring scale as the one used to evaluate the frost damages in the field experiment, except that scores was attributed to single plants instead of plots.

Overall, 5 traits constituted the phenotyping data for the GWAS, abbreviated as follows: FD_CC: frost damages in the controlled conditions experiment, FD_Field_date1, FD_Field_date2, FD_Field_date3 and FD_Field_date4: winter frost damages evaluated in the field experiment at the first, second, third and fourth date respectively.

## Phenotypic data analyses

The phenotypic data were analyzed with the R 3.5.0 software [56, 57] using a linear mixed model to obtain estimates of variance components, heritability ($H^2$), as well as best linear unbiased predictions (BLUPs) of adjusted means. The following Linear Mixed Model (LMM) was used: $Y_{ij} = \mu + geno_i + rep_j + e_{ij}$, where $Y_{ij}$ is the value of frost damages recorded for the genotype *i* at the replicate *j*. $\mu$ is the mean, $geno_i$ is the random genetic effect of the genotype *i*, $rep_j$ is the fixed replicate effect of the replicate *j* and $e_{ij}$ is the residual effect. The model was carried out using the R function “lmer” of the package ‘lme4’ [58, 59]. Heritability ($H^2$) was estimated using the following formula: $H^2 = V_g/(V_g + V_{res}/n_{rep})$, where $V_g$ is

the genotypic variance component, $V_{res}$ is the residual variance component and $n_{rep}$ is the number of replicates taking into account the missing values. BLUPs for each genotype-trait combination were calculated from each LMM analysis using the function “ranef”, implemented in ‘lme4’ package of R software [59] and were used for the GWA analysis.

## Genotyping and quality control

The collection was genotyped at 11,366 SNPs using the Illumina Infinium® BeadChip 13.2 K SNPs as described in [36]. These SNPs were all located in gene-context sequences and derived from separated transcripts [60]. The consensus genetic map from Tayeh et al. [36] was used as the genetic framework for the association analyses. This map was built on the basis of genotyping data collected for 12 pea recombinant inbred line populations. Considering the large collinearity between individual maps, a set of genotyping data for 15,352 markers

and from all populations was used to build the consensus map. The latter shows a cumulative total length of 794.9 cM and a mean inter-marker distance of 0.24 cM.

The genotyping matrix, which was composed of a set of 11,366 SNPs and 365 pea accessions, was filtered using Plink v.1.9 software [61, 62]. Accessions and SNP markers with a call rate below 0.90 as well as SNP markers with a minor allele frequency (MAF) below 0.05 were excluded from the GWA analysis. After quality control checking, a genotyping matrix consisting of 10, 739 SNPs and 363 accessions with 0.6% missing data was kept for further analyses. The resulting data set was further imputed using Beagle v.3.3.2 software [63]. Beagle applies a Markov model to the hidden states (the haplotype phase and the true genotype) along the chromosome using an EM (Expectation-Maximization) algorithm that iteratively updates model parameters to maximize the model likelihood up to the moment where convergence is achieved. Finally, a genotyping matrix consisting of 10,739 SNPs and 363 accessions with no missing data was used for GWAS. Scripts from Negro et al. [64] were used for the quality control and imputation.

## Linkage disequilibrium estimation

The estimates of linkage disequilibrium (LD) within the collection were determined by the squared allelic-frequency correlations ($r^2$) for pairs of loci as described in Weir [65]. Linkage disequilibrium analysis between pairs of SNP markers was calculated in a sliding window of 900 markers using Plink v.1.9 software [61, 62]. Then, intrachromosomal LD quantification and graphical representation of LD decay were accomplished using R 3.5.0 software [56, 57]. The LD decay was measured as the genetic distance (cM) where the average $r^2$ decreased to half its maximum value.

## Population structure and individual relatedness

To control false positive associations, population structure and individual relatedness (kinship) among accessions of the collection were taken into account by fitting markers based structure and kinship matrices in the association models [66]. Kinship and population structure were estimated using a matrix data composed of 363 accessions and a set of 2962 markers without any missing data and corresponding to non-redundant genetic positions randomly selected on the consensus map. The coefficients of kinship between pairs of accessions were estimated using the realized relationship matrix kinship estimation approach implemented in FaST-LMM software [67]. Two alternative approaches were considered to estimate the kinship matrix as described by Rincent et al. [68]. In the first one, the kinship was estimated with all the markers that are not located on the same linkage group (LG) than the tested SNP. Thus,
environment chamber. The temperature, light level and humidity were recorded and were similar during the three experiments. Briefly, the plants at the stage of 2nd - 3rd leaf were first treated with a regime of 11 days of cold acclimation at 10 °C/2 °C (day/night) with a 10 h photoperiod. The frost treatment was then carried out at 6 °C/-8 °C with 8 h of daylight during 4 days. After frost, a recovery period was applied with a temperature regime of 16 °C/5 °C and 10 h of daylight during 8 days. Frost tolerance was evaluated by scoring frost damages at the end of the recovery period with the same scoring scale as the one used to evaluate the frost damages in the field experiment, except that scores was attributed to single plants instead of plots.

Overall, 5 traits constituted the phenotyping data for the GWAS, abbreviated as follows: FD_CC: frost damages in the controlled conditions experiment, FD_Field_date1, FD_Field_date2, FD_Field_date3 and FD_Field_date4: winter frost damages evaluated in the field experiment at the first, second, third and fourth date respectively.

## Phenotypic data analyses

The phenotypic data were analyzed with the R 3.5.0 software [56, 57] using a linear mixed model to obtain estimates of variance components, heritability ($H^2$), as well as best linear unbiased predictions (BLUPs) of adjusted means. The following Linear Mixed Model (LMM) was used: $Y_{ij} = \mu + geno_i + rep_j + e_{ij}$, where $Y_{ij}$ is the value of frost damages recorded for the genotype *i* at the replicate *j*. $\mu$ is the mean, $geno_i$ is the random genetic effect of the genotype *i*, $rep_j$ is the fixed replicate effect of the replicate *j* and $e_{ij}$ is the residual effect. The model was carried out using the R function “lmer” of the package ‘lme4’ [58, 59]. Heritability ($H^2$) was estimated using the following formula: $H^2 = V_g/(V_g + V_{res}/n_{rep})$, where $V_g$ is the genotypic variance component, $V_{res}$ is the residual variance component and $n_{rep}$ is the number of replicates taking into account the missing values. BLUPs for each genotype-trait combination were calculated from each LMM analysis using the function “ranef”, implemented in ‘lme4’ package of R software [59] and were used for the GWA analysis.

## Genotyping and quality control

The collection was genotyped at 11,366 SNPs using the Illumina Infinium® BeadChip 13.2 K SNPs as described in [36]. These SNPs were all located in gene-context sequences and derived from separated transcripts [60]. The consensus genetic map from Tayeh et al. [36] was used as the genetic framework for the association analyses. This map was built on the basis of genotyping data collected for 12 pea recombinant inbred line populations. Considering the large collinearity between individual maps, a set of genotyping data for 15,352 markers

and from all populations was used to build the consensus map. The latter shows a cumulative total length of 794.9 cM and a mean inter-marker distance of 0.24 cM.

The genotyping matrix, which was composed of a set of 11,366 SNPs and 365 pea accessions, was filtered using Plink v.1.9 software [61, 62]. Accessions and SNP markers with a call rate below 0.90 as well as SNP markers with a minor allele frequency (MAF) below 0.05 were excluded from the GWA analysis. After quality control checking, a genotyping matrix consisting of 10, 739 SNPs and 363 accessions with 0.6% missing data was kept for further analyses. The resulting data set was further imputed using Beagle v.3.3.2 software [63]. Beagle applies a Markov model to the hidden states (the haplotype phase and the true genotype) along the chromosome using an EM (Expectation-Maximization) algorithm that iteratively updates model parameters to maximize the model likelihood up to the moment where convergence is achieved. Finally, a genotyping matrix consisting of 10,739 SNPs and 363 accessions with no missing data was used for GWAS. Scripts from Negro et al. [64] were used for the quality control and imputation.

## Linkage disequilibrium estimation

The estimates of linkage disequilibrium (LD) within the collection were determined by the squared allelic-frequency correlations ($r^2$) for pairs of loci as described in Weir [65]. Linkage disequilibrium analysis between pairs of SNP markers was calculated in a sliding window of 900 markers using Plink v.1.9 software [61, 62]. Then, intrachromosomal LD quantification and graphical representation of LD decay were accomplished using R 3.5.0 software [56, 57]. The LD decay was measured as the genetic distance (cM) where the average $r^2$ decreased to half its maximum value.

## Population structure and individual relatedness

To control false positive associations, population structure and individual relatedness (kinship) among accessions of the collection were taken into account by fitting markers based structure and kinship matrices in the association models [66]. Kinship and population structure were estimated using a matrix data composed of 363 accessions and a set of 2962 markers without any missing data and corresponding to non-redundant genetic positions randomly selected on the consensus map. The coefficients of kinship between pairs of accessions were estimated using the realized relationship matrix kinship estimation approach implemented in FaST-LMM software [67]. Two alternative approaches were considered to estimate the kinship matrix as described by Rincent et al. [68]. In the first one, the kinship was estimated with all the markers that are not located on the same linkage group (LG) than the tested SNP. Thus,
seven kinship matrices were estimated, each being specific to a linkage group; these matrices were noted K<sub>LGx</sub> with *x* corresponding to the number of linkage group tested. Such an approach aims at increasing power of detection of significant markers in GWAS particularly in regions of high LD. In the second approach, correlation between markers took into account all the 2962 markers and the kinship matrix was noted K.

The discriminant analysis of principal components (DAPC) method developed by Jombard et al. [69] and implemented into the ‘adegenet’ R package [70–72] was used to cluster accessions on the basis of their genotype. This method aims at identifying and describing clusters of genetically related individuals without prior knowledge of groups. First, the optimal number of genetic clusters (k) was determined through the ‘K-means’ method using the function “find.clusters”. The number of clusters was allowed to vary from one to 20 during the determination of the optimal value of k, based on the Bayesian Information Criterion (BIC). The most likely number of clusters was chosen on the basis of the lowest associated BIC. Then, the principal component analysis (PCA) step of DAPC was performed through maximization of the ‘a-score’ criterion and the optimal number of principal components (PCs) was obtained after 100 iterations using the function “optim.a.score” implemented in ‘adegenet’ package of R software. Finally, DAPC was performed considering the most likely number of clusters (k) and the optimal number of PCs identified using the function “dapc” implemented in adegenet R-package [70–72]. To confirm the allocation of accessions to clusters by DAPC analysis, a Nei genetic distance matrix [73] was calculated with the function “stampp-NeisD” implemented in ‘StAMPP’ package of R software [74] using the genotyping data composed of 363 accessions and a set of 2962 SNPs. Then the resulting matrix was plotted as a dendrogram using the ward method with the package ‘cluster’ implemented in R software [75].

The two first coordinates of DAPC results were used as covariates (Q matrix) in the GWAS to correct the association tests for false positives.

## Association mapping

BLUPs corresponding to the phenotypic data collected for each accession were used to identify marker-trait association using Linear Mixed Model (LMM) accounting for kinship matrix (K or K<sub>LG</sub>) with or without population structure matrix (Q) as random effect(s). Four models were therefore compared for their capacity to fit the data: (1) a LMM using the kinship matrix K, (2) a LMM corrected for kinship using the K<sub>LG</sub> matrices, (3) a LMM including the K and Q matrices and (4) a LMM using both K<sub>LG</sub> and Q matrices. For each frost damages trait, the best model was chosen by comparing the

likelihoods of each model using the Bayesian Information Criterion (BIC) [76]. The model with the smallest BIC was selected. All analyses were performed using LMM provided by the FaST-LMM version algorithm [67]. The threshold to declare an association significant was set at a probability level of the *p*-value inferior to 4.65E-06, i.e. -log10 (*p*) > 5.33, which corresponded to the Bonferroni threshold (0.05/ number of tested SNPs). To represent the association results, Manhattan plots and their corresponding quantile-quantile plots were drawn using the package ‘QQman’ implemented in R software [77].

## Local LD block estimation and favorable allele identification

Local LD analysis was used to define the LD blocks around significant associated markers detected by GWAS using Plink 1.9 software [61, 62]. For each associated marker, markers in strong linkage disequilibrium (LD; *r*² > 0.8) with this one, were identified to define a LD block. By this way, a LD block was defined as the interval including all markers in LD (*r*² > 0.8) with the targeted associated marker(s). Unique associated markers which didn’t constitute a LD block were kept for further analyses. Thus, for each identified genomic region, LD blocks and unique associated marker(s) composed a significant locus related with frost tolerance. At each significant locus, haplotypes were identified, among the accessions of the collection, according to the non-imputed genotyping data corresponding to the list of markers significantly detected by GWAS and linked markers. Haplotypes showing missing data loci as well as SNP with heterozygous genotypic data were excluded from further analysis. Besides, haplotypes represented by less than 5% of the total number of accessions were also removed from the analysis. Based on the results of association mapping, the allelic effect corresponding to the minor allele (a<sub>eff</sub>) of markers significantly associated with the frost damage traits were analyzed: if a<sub>eff</sub> had a negative value, the minor allele of the associated marker was considered to decrease frost damage (favorable allele for frost tolerance); if a<sub>eff</sub> had a positive value, the minor allele of the associated marker was considered to increase frost damage (unfavorable allele for frost tolerance). For each significant locus and each corresponding trait, the values of frost damages of the different haplotypes were compared using an analysis of variance with a nested design for ‘haplotype/genotype’, followed by a Student-Newman-Keuls (SNK) comparison test using the function “SNK.test” of the R-package ‘agricolae’ [78]. Favorable and unfavorable haplotypes at each significant locus were defined as follows: the favorable haplotypes should show a significant lower frost damage mean score while unfavorable haplotypes should show a significant higher frost
seven kinship matrices were estimated, each being specific to a linkage group; these matrices were noted K<sub>LGx</sub> with *x* corresponding to the number of linkage group tested. Such an approach aims at increasing power of detection of significant markers in GWAS particularly in regions of high LD. In the second approach, correlation between markers took into account all the 2962 markers and the kinship matrix was noted K.

The discriminant analysis of principal components (DAPC) method developed by Jombard et al. [69] and implemented into the ‘adegenet’ R package [70–72] was used to cluster accessions on the basis of their genotype. This method aims at identifying and describing clusters of genetically related individuals without prior knowledge of groups. First, the optimal number of genetic clusters (k) was determined through the ‘K-means’ method using the function “find.clusters”. The number of clusters was allowed to vary from one to 20 during the determination of the optimal value of k, based on the Bayesian Information Criterion (BIC). The most likely number of clusters was chosen on the basis of the lowest associated BIC. Then, the principal component analysis (PCA) step of DAPC was performed through maximization of the ‘a-score’ criterion and the optimal number of principal components (PCs) was obtained after 100 iterations using the function “optim.a.score” implemented in ‘adegenet’ package of R software. Finally, DAPC was performed considering the most likely number of clusters (k) and the optimal number of PCs identified using the function “dapc” implemented in adegenet R-package [70–72]. To confirm the allocation of accessions to clusters by DAPC analysis, a Nei genetic distance matrix [73] was calculated with the function “stampp-NeisD” implemented in ‘StAMPP’ package of R software [74] using the genotyping data composed of 363 accessions and a set of 2962 SNPs. Then the resulting matrix was plotted as a dendrogram using the ward method with the package ‘cluster’ implemented in R software [75].

The two first coordinates of DAPC results were used as covariates (Q matrix) in the GWAS to correct the association tests for false positives.

## Association mapping

BLUPs corresponding to the phenotypic data collected for each accession were used to identify marker-trait association using Linear Mixed Model (LMM) accounting for kinship matrix (K or K<sub>LG</sub>) with or without population structure matrix (Q) as random effect(s). Four models were therefore compared for their capacity to fit the data: (1) a LMM using the kinship matrix K, (2) a LMM corrected for kinship using the K<sub>LG</sub> matrices, (3) a LMM including the K and Q matrices and (4) a LMM using both K<sub>LG</sub> and Q matrices. For each frost damages trait, the best model was chosen by comparing the

likelihoods of each model using the Bayesian Information Criterion (BIC) [76]. The model with the smallest BIC was selected. All analyses were performed using LMM provided by the FaST-LMM version algorithm [67]. The threshold to declare an association significant was set at a probability level of the *p*-value inferior to 4.65E-06, i.e. -log10 (*p*) > 5.33, which corresponded to the Bonferroni threshold (0.05/ number of tested SNPs). To represent the association results, Manhattan plots and their corresponding quantile-quantile plots were drawn using the package ‘QQman’ implemented in R software [77].

## Local LD block estimation and favorable allele identification

Local LD analysis was used to define the LD blocks around significant associated markers detected by GWAS using Plink 1.9 software [61, 62]. For each associated marker, markers in strong linkage disequilibrium (LD; *r*² > 0.8) with this one, were identified to define a LD block. By this way, a LD block was defined as the interval including all markers in LD (*r*² > 0.8) with the targeted associated marker(s). Unique associated markers which didn’t constitute a LD block were kept for further analyses. Thus, for each identified genomic region, LD blocks and unique associated marker(s) composed a significant locus related with frost tolerance. At each significant locus, haplotypes were identified, among the accessions of the collection, according to the non-imputed genotyping data corresponding to the list of markers significantly detected by GWAS and linked markers. Haplotypes showing missing data loci as well as SNP with heterozygous genotypic data were excluded from further analysis. Besides, haplotypes represented by less than 5% of the total number of accessions were also removed from the analysis. Based on the results of association mapping, the allelic effect corresponding to the minor allele (a<sub>eff</sub>) of markers significantly associated with the frost damage traits were analyzed: if a<sub>eff</sub> had a negative value, the minor allele of the associated marker was considered to decrease frost damage (favorable allele for frost tolerance); if a<sub>eff</sub> had a positive value, the minor allele of the associated marker was considered to increase frost damage (unfavorable allele for frost tolerance). For each significant locus and each corresponding trait, the values of frost damages of the different haplotypes were compared using an analysis of variance with a nested design for ‘haplotype/genotype’, followed by a Student-Newman-Keuls (SNK) comparison test using the function “SNK.test” of the R-package ‘agricolae’ [78]. Favorable and unfavorable haplotypes at each significant locus were defined as follows: the favorable haplotypes should show a significant lower frost damage mean score while unfavorable haplotypes should show a significant higher frost
damage mean score. Finally, we listed representative accessions for each favorable and unfavorable haplotype based on the following condition: each accession should show a mean score of the considered associated trait(s) inferior to 1 for favorable haplotypes and superior to 4 for unfavorable ones.

## Annotated genes underlying frost tolerance loci

To identify genes that may be associated with the frost damage phenotypes, a region encompassing 1 Mb flanking regions upstream and downstream from each of the FD-related markers, ie. significant GWAS markers and markers in LD ($r^2$ > 0.8) with the former ones, was defined. This region was searched for genes annotated in the pea genome assembly v.1a developed by Kreplak et al. [39] using the genome JBrowse available at https://urgi.versailles.inra.fr/Species/Pisum.

## Supplementary information

Supplementary information accompanies this paper at https://doi.org/10.1186/s12864-020-06928-w.

**Additional file 1: Figure S1.** Distribution of Best Linear Unbiased Prediction (BLUP) values for the five traits observed within the pea collection. A: frost damages in the controlled conditions experiment. B, C, D and E: frost damages in the field experiment at the date 1, 2, 3 and 4 respectively.

**Additional file 2: Figure S2.** Distribution of 10,739 SNPs along the *Pisum sativum* linkage groups. Number of SNPs per position are indicated as grey horizontal bars. Genetic position in cM is shown on the y-axis and number of SNPs per position is shown on the x-axis.

**Additional file 3: Table S1.** Description of the *Pisum sativum* linkage groups (LGs) used in the present study. The number of SNP markers, the genetic length (in cM, from Tayeh et al. [36]) and the average minor allele frequency (MAF) are shown for each LG.

**Additional file 4: Figure S3.** Distribution of minor allele frequencies (MAF) for 10,739 SNP markers within the 363 pea accessions.

**Additional file 5: Figure S4.** Scatterplot showing the linkage disequilibrium (LD) decay estimated in the association mapping collection. The LD decay across each linkage group (LG) and the overall LD decay across the genome (All LG) are shown. The $r^2$ values of LD between pairs of markers considered are plotted as a function of the genetic position in cM. Red curves represent the estimated LD decay. Blue dashed horizontal lines represent half of the maximum LD value. Blue dashed vertical lines represent the estimated genetic distance (cM) at which the LD decay dropped to half of its maximum. LD decay rate is represented as the point of intersection between the two dashed lines.

**Additional file 6: Table S2.** Description of the association mapping collection. This table presents the list of the pea accessions composing the association mapping collection with their end-use, cultivation status, geographical origin and sowing type. The ‘DAPC Cl Cluster (k)’ column shows assignation of the 363 pea accessions to a cluster based on the discriminant analysis of principal components (DAPC). The ‘Dendrogram Cl Cluster (k)’ column shows the allocation of individuals to clusters based on the dendrogram using Nei genetic distances between accessions. The description of the pea accessions is extracted from Burstin et al. [55]. CRB Code: Code used for the association mapping, also named collection of biological resources (CRB). *: sowing type modified for this accession, according to Yenne et al. [79].

**Additional file 7: Figure S5.** Dendrogram from Nei genetic distance matrix for 363 genotypes of the pea reference collection. On the y-axis are represented the genetic distances between clusters or accessions. On

the x-axis are represented, in red font, the clusters identified for a Nei genetic distance of 7.

**Additional file 8: Figure S6.** Distribution of the kinship coefficients between accessions of the association mapping collection. The first histogram (A) describes the distribution of the kinship coefficients within the K matrix, calculated with all markers of the genome. The remaining histograms (B, C, D, E, F, G and H) describe the kinship coefficients within each of the seven $K_{LG}$ matrices calculated as explained in the material and methods section (for example the kinship matrix $K_{LG1}$ was estimated with all the markers except those that are located on the first linkage group).

**Additional file 9: Table S3.** Description of linkage disequilibrium (LD) blocks per linkage group in the association mapping collection. A LD block consists in a series of at least 2 markers which are in significant LD ($r^2$ > 0.8) with at least one trait-associated marker (underlined marker). LD blocks are named in consecutive numerical order following their linkage group (LG) name. cM*: genetic position, in centiMorgan of each marker along the genetic map of the corresponding linkage group; LD ($r^2$) **: $r^2$ value of each marker with the other markers of the same LD block.

**Additional file 10: Table S4.** Marker haplotype analysis of the association mapping collection. For each linkage group (LG), the list of markers significantly detected by GWAS and markers in linkage disequilibrium (LD); $r^2$ > 0.8) with the former ones is shown. The third line shows genetic positions from the consensus map of Tayeh et al. [36]. The fourth line indicates the LD blocks composed and named as mentioned in the legend of Additional file 9. The following lines show the allelic composition of haplotypes defined by LD blocks and individual associated markers at each of the 6 frost tolerance loci on linkage groups (LGs) I, II, III, V, VI and VII. For each frost damage (FD)-associated marker, the favorable allele is in red font and the unfavorable allele in blue font. Haplotypes are named in consecutive numerical order following their linkage group name; only haplotypes without missing values or heterozygous markers and carried by more than 3% of the lines from the association mapping collection are listed. For each haplotype, accessions and their mean phenotypic values ± standard error of the variables significantly associated with marker(s) in the linkage group are shown. Significant differences between haplotypes were assayed by a SNK means comparison test; favorable haplotypes are shown by a red background and unfavorable haplotypes by a blue background, regarding the SNK test. Haplotypes with a white background are classified in intermediate groups.

**Additional file 11: Table S5.** List of annotated genes underlying genome-wide association loci of frost tolerance in pea. Genes that were located in an interval of ±1 Mb on both sides of markers significantly detected by GWAS and markers in linkage disequilibrium (LD); $r^2$ > 0.8) with the former ones, are listed. For each identified gene, the nearest marker significantly detected by GWAS (underlined font) or marker in LD with associated marker(s) (non-underlined font) is shown. The annotation of genes was extracted from *Pisum sativum* v.1a genome JBrowse available at https://urgi.versailles.inra.fr/Species/Pisum [39]. Genes positions in the pea genome assembly v.1a are presented by their assigned chromosomes and physical positions indicated in bp. Genes which are not assigned to one of the seven chromosomes, are represented by their physical positions on unanchored scaffolds.ª: Annotation refined with the homologous gene from *Medicago truncatula* available on the *Pisum sativum* v.1a genome JBrowse, and whose corresponding gene function was identified from the *Medicago truncatula* v4.0 genome JBrowse available on www.medicagogenome.org.ᵇ: Annotation refined with the homologous gene from *Glycine max* available on the *Pisum sativum* v.1a genome JBrowse and whose corresponding gene function was identified from the genome v9.0 assembly V1.1 available on http://soykb.org/gene supra.php.ᶜ: Annotation refined with the predicted protein function of transcript sequences corresponding to mapped SNPs, extracted from Table S10 in Tayeh et al. [36].

**Additional file 12: Table S6.** Description of the pea association mapping collection as described in Additional file 6 and the correspondence of the genotyping results of accessions at the *Hr* locus.
damage mean score. Finally, we listed representative accessions for each favorable and unfavorable haplotype based on the following condition: each accession should show a mean score of the considered associated trait(s) inferior to 1 for favorable haplotypes and superior to 4 for unfavorable ones.

## Annotated genes underlying frost tolerance loci

To identify genes that may be associated with the frost damage phenotypes, a region encompassing 1 Mb flanking regions upstream and downstream from each of the FD-related markers, ie. significant GWAS markers and markers in LD ($r^2$ > 0.8) with the former ones, was defined. This region was searched for genes annotated in the pea genome assembly v.1a developed by Kreplak et al. [39] using the genome JBrowse available at https://urgi.versailles.inra.fr/Species/Pisum.

## Supplementary information

Supplementary information accompanies this paper at https://doi.org/10.1186/s12864-020-06928-w.

**Additional file 1: Figure S1.** Distribution of Best Linear Unbiased Prediction (BLUP) values for the five traits observed within the pea collection. A: frost damages in the controlled conditions experiment. B, C, D and E: frost damages in the field experiment at the date 1, 2, 3 and 4 respectively.

**Additional file 2: Figure S2.** Distribution of 10,739 SNPs along the *Pisum sativum* linkage groups. Number of SNPs per position are indicated as grey horizontal bars. Genetic position in cM is shown on the y-axis and number of SNPs per position is shown on the x-axis.

**Additional file 3: Table S1.** Description of the *Pisum sativum* linkage groups (LGs) used in the present study. The number of SNP markers, the genetic length (in cM, from Tayeh et al. [36]) and the average minor allele frequency (MAF) are shown for each LG.

**Additional file 4: Figure S3.** Distribution of minor allele frequencies (MAF) for 10,739 SNP markers within the 363 pea accessions.

**Additional file 5: Figure S4.** Scatterplot showing the linkage disequilibrium (LD) decay estimated in the association mapping collection. The LD decay across each linkage group (LG) and the overall LD decay across the genome (All LG) are shown. The $r^2$ values of LD between pairs of markers considered are plotted as a function of the genetic position in cM. Red curves represent the estimated LD decay. Blue dashed horizontal lines represent half of the maximum LD value. Blue dashed vertical lines represent the estimated genetic distance (cM) at which the LD decay dropped to half of its maximum. LD decay rate is represented as the point of intersection between the two dashed lines.

**Additional file 6: Table S2.** Description of the association mapping collection. This table presents the list of the pea accessions composing the association mapping collection with their end-use, cultivation status, geographical origin and sowing type. The ‘DAPC_Cluster (k)’ column shows assignation of the 363 pea accessions to a cluster based on the discriminant analysis of principal components (DAPC). The ‘Dendrogram_Cluster (k)’ column shows the allocation of individuals to clusters based on the dendrogram using Nei genetic distances between accessions. The description of the pea accessions is extracted from Burstin et al. [55]. CRB Code: Code used for the association mapping, also named collection of biological resources (CRB). *: sowing type modified for this accession, according to Yenne et al. [79].

**Additional file 7: Figure S5.** Dendrogram from Nei genetic distance matrix for 363 genotypes of the pea reference collection. On the y-axis are represented the genetic distances between clusters or accessions. On

the x-axis are represented, in red font, the clusters identified for a Nei genetic distance of 7.

**Additional file 8: Figure S6.** Distribution of the kinship coefficients between accessions of the association mapping collection. The first histogram (A) describes the distribution of the kinship coefficients within the K matrix, calculated with all markers of the genome. The remaining histograms (B, C, D, E, F, G and H) describe the kinship coefficients within each of the seven $K_{LG}$ matrices calculated as explained in the material and methods section (for example the kinship matrix $K_{LG1}$ was estimated with all the markers except those that are located on the first linkage group).

**Additional file 9: Table S3.** Description of linkage disequilibrium (LD) blocks per linkage group in the association mapping collection. A LD block consists in a series of at least 2 markers which are in significant LD ($r^2$ > 0.8) with at least one trait-associated marker (underlined marker). LD blocks are named in consecutive numerical order following their linkage group (LG) name. cM*: genetic position, in centiMorgan of each marker along the genetic map of the corresponding linkage group; LD ($r^2$) **: $r^2$ value of each marker with the other markers of the same LD block.

**Additional file 10: Table S4.** Marker haplotype analysis of the association mapping collection. For each linkage group (LG), the list of markers significantly detected by GWAS and markers in linkage disequilibrium (LD); $r^2$ > 0.8) with the former ones is shown. The third line shows genetic positions from the consensus map of Tayeh et al. [36]. The fourth line indicates the LD blocks composed and named as mentioned in the legend of Additional file 9. The following lines show the allelic composition of haplotypes defined by LD blocks and individual associated markers at each of the 6 frost tolerance loci on linkage groups (LGs) I, II, III, V, VI and VII. For each frost damage (FD)-associated marker, the favorable allele is in red font and the unfavorable allele in blue font. Haplotypes are named in consecutive numerical order following their linkage group name; only haplotypes without missing values or heterozygous markers and carried by more than 3% of the lines from the association mapping collection are listed. For each haplotype, accessions and their mean phenotypic values ± standard error of the variables significantly associated with marker(s) in the linkage group are shown. Significant differences between haplotypes were assayed by a SNK means comparison test; favorable haplotypes are shown by a red background and unfavorable haplotypes by a blue background, regarding the SNK test. Haplotypes with a white background are classified in intermediate groups.

**Additional file 11: Table S5.** List of annotated genes underlying genome-wide association loci of frost tolerance in pea. Genes that were located in an interval of ±1 Mb on both sides of markers significantly detected by GWAS and markers in linkage disequilibrium (LD); $r^2$ > 0.8) with the former ones, are listed. For each identified gene, the nearest marker significantly detected by GWAS (underlined font) or marker in LD with associated marker(s) (non-underlined font) is shown. The annotation of genes was extracted from *Pisum sativum* v.1a genome JBrowse available at https://urgi.versailles.inra.fr/Species/Pisum [39]. Genes positions in the pea genome assembly v.1a are presented by their assigned chromosomes and physical positions indicated in bp. Genes which are not assigned to one of the seven chromosomes, are represented by their physical positions on unanchored scaffolds.ª: Annotation refined with the homologous gene from *Medicago truncatula* available on the *Pisum sativum* v.1a genome JBrowse, and whose corresponding gene function was identified from the *Medicago truncatula* v4.0 genome JBrowse available on www.medicagogenome.org.ᵇ: Annotation refined with the homologous gene from *Glycine max* available on the *Pisum sativum* v.1a genome JBrowse and whose corresponding gene function was identified from the genome v9.0 assembly V1.1 available on http://soykb.org/gene supra.php.ᶜ: Annotation refined with the predicted protein function of transcript sequences corresponding to mapped SNPs, extracted from Table S10 in Tayeh et al. [36].

**Additional file 12: Table S6.** Description of the pea association mapping collection as described in Additional file 6 and the correspondence of the genotyping results of accessions at the *Hr* locus.
**Additional file 13: Table S7.** Description of accessions of the pea reference collection for their haplotype at the frost damage (FD)-associated loci of the GWA study and at the *Hr* locus. At the linkage groups (LGs) I, II, III, V, VI and VII, the favorable haplotypes are shown by a red background, the unfavorable haplotypes by a blue background, as described in the legend of Additional file 10. Accessions with undefined haplotypes or intermediate haplotypes were not presented. The same colour code has been used to describe the favorable (red background: *Hr*) and unfavorable (blue background: *hr*) allele for the *Hr* gene, as determined in Lejeune-Hénaut et al. [21]. The mean frost damage score observed in the field experiment as well as its standard error (SE) are given. Frost scores are ranging from 0 (no damage) to 5 (dead plant). The passport data of accessions are extracted from the Additional file 6.

**Additional file 14: Table S8.** Daily air and soil temperatures during the field experiment in Clermont-Ferrand Theix.

## Abbreviations

AP2/EREBP: APETALA2/Ethylene Responsive Element Binding Protein;
BIC: Bayesian Information Criterion; BLUPs: Best Linear Unbiased Predictions;
CBF: C-repeat Binding Factor; cM: centi Morgan; COR: Cold Regulated;
CVg: Coefficient of Variation of the genetic variance; DAPC: Discriminant
Analysis of Principal Components; DREB: Dehydration-Responsive Element-
Binding; Elf3: Early flowering 3; FD: Frost Damage; FD_CC: Frost damages in
the controlled conditions experiment; FD_Field_date1: Frost damages in the
field experiment at the first date of scoring; FD_Field_date2: Frost damages
in the field experiment at the second date of scoring; FD_Field_date3: Frost
damages in the field experiment at the third date of scoring;
FD_Field_date4: Frost damages in the field experiment at the fourth date of
scoring; GWAS: Genome Wide Association Study; H²: Broad sense heritability;
Hr: High response to photoperiod; INRAE: National Research Institute for
Agriculture, Food and Environment; k: DAPC Cluster; K: Kinship matrix;
K_LG: Kinship matrix corresponding to a given Linkage Group; LD: Linkage
Disequilibrium; LG: Linkage Group; LMM: Linear Mixed Model; MAF: Minor
Allele Frequency; NAM: Nested Association Mapping; NIL: Near Isogenic Line;
PCs: Principal Components; PCA: Principal Components Analysis; Q: Structure
matrix; QTL: Quantitative Trait Locus; RIL: Recombinant Inbred Line;
SE: Standard Error; SNK: Student-Newman-Keuls; SNP: Single Nucleotide
Polymorphism; SSD: Single Seed Descent; SSR: Simple Sequence Repeat;
Vg: Genetic variance; WFD: Winter Frost Damage

## Acknowledgements

The authors would like to thank the INRAE experimental units of Estrées-Mons, Dijon and Theix, France, for their contribution to the field experiments. They are also grateful to Catherine Desmetz (Agroécologie, INRAE) for additional genotyping of the reference collection, to Marianne Chabert-Martinello, Anthony Klein and Jean-Bernard Magnin-Robert (Agroécologie, INRAE) for providing informations about the constitution of the reference collection and to Renaud Rincent (GDEC, INRAE) for his advices in genome wide association analyses.

## Authors' contributions

VF, RD and MT collected phenotypic data in the field and controlled conditions experiments. GA and JB coordinated the genotyping of the association mapping collection and provided genotyping data. SB performed the phenotypic and genetic analyses, as well as association mapping and haplotypes analysis and wrote the manuscript. MS and SSN assisted in analysing genetic data by testing GWAS models and providing a part of R scripts, respectively. ILH and BD conceived and coordinated the study. ILH managed the technical and financial supplies of the study within the PeaMUST project. JLH, BD and ILH coordinated the overall progress and financial support of the study. ILH, BD and NB reviewed and contributed to draft the manuscript. All authors read and approved the final manuscript.

## Funding

This work was supported by the PhD fellowship of Sana Beji which was co-funded by the University of Lille, doctoral school for materials, radiation and environmental sciences (France) and by the region Hauts-de-France (France). The acquisition of the phenotyping data was funded by the SNPEA project (ANR06-GPLA019). The acquisition of the genotyping data was funded by the GENOPEA project (ANR-09-GENM-026). The analyzes, interpretation of

data and writing of the manuscript were realized in the frame of the Pea-MUST project (ANR11-BTBR-0002). The three project fundings come from the French government and were managed by the Research National Agency (ANR). The funding agency did not contribute to the design of the study or collection, nor to the analysis and interpretation of data or to the writing of the manuscript.

## Availability of data and materials

Genetic positions and sequences of markers used for the current study are available in Tayeh et al. [36]. The pea genome assembly v.1a explored in this study is available at https://urgi.versailles.inra.fr/Species/Pisum (Kreplak et al. [39]). The data supporting the findings of this study were used under license within the PeaMUST project and are not publicly available. Phenotyping and genotyping data, all other intermediate data files and scripts are however available from the authors upon reasonable request and subjected to data transfer agreement.

## Ethics approval and consent to participate

Not applicable.

## Consent for publication

Not applicable.

## Competing interests

The authors declare that they have no competing interests.

## Author details

¹BioEcoAgro, INRAE, Univ. Liège, Univ. Lille, Univ. Picardie Jules Verne, 2, Chaussée Brunehaut, F-80203 Estrées-Mons, France. ²GCIE-Piccardie, INRAE, F-80203 Estrées-Mons, France. ³GQE - Le Moulon, INRAE, Univ. Paris-Sud, CNRS, AgroParisTech, Univ. Paris-Saclay, F-91190 Gif-sur-Yvette, France.

⁴Agroécologie, AgroSup Dijon, INRAE, Univ. Bourgogne, Univ. Bourgogne Franche-Comté, F-21000 Dijon, France.

Received: 9 August 2019 Accepted: 20 July 2020

Published online: 04 August 2020

## References

1. FAOSTAT. FAOSTAT data base. Food and agriculture organization of the United Nations - Statistic Division. 2018. http://www.fao.org/faostat/en/#data/. Accessed 27 Mar 2020.

2. Graham P, Vance C. Legumes: importance and constraints to greater use. Plant Physiol. 2003;131:872-7. https://doi.org/10.1104/pp.017004.

3. Jeuffroy MH, Baranger E, Carrouee B, de Chezelles E, Gosme M, Henault C, et al. Nitrous oxide emissions from crop rotations including wheat, oilseed rape and dry peas. Biogeosciences. 2013;10(3):1787-97. https://doi.org/10.5194/bg-10-1787-2013.

4. Benezit M, Biarnes V, Jeuffroy MH. Impact of climate and diseases on pea yields: what perspectives with climate change? OCL. 2017;24(1):D103. https://doi.org/10.1051/ocl/2016055.

5. Lejeune-Hénaut I, Bourion V, Eteve G, Cunot E, Delhaye K, Desmyter C. Floral initiation in field-grown forage peas is delayed to a greater extent by short photoperiods, than in other types of European varieties. Euphytica. 1999; 109(3):201-11. https://doi.org/10.1023/a:1003727324475.

6. Fowler DB, Breton G, Limin AE, Mahfoozi S, Sarhan F. Photoperiod and temperature interactions regulate low-temperature-induced gene expression in barley. Plant Physiol. 2001;127(4):1676-81. https://doi.org/10.1104/pp.010483.

7. Thomashow M. Plant cold acclimation, freezing tolerance genes and regulatory mechanisms. Annu Rev Plant Physiol Plant Mol Biol. 1999;50:571-99. https://doi.org/10.1146/annurev.arplant.50.1.571.

8. Ruelland E, Vaultier M, Zachowski A, Hurry V. Cold signalling and cold acclimation in plants. Adv Bot Res. 2009;49:35-146. https://doi.org/10.1016/S0065-2296(08)00602-2.

9. Rihan HZ, Al-Issawi M, Fuller MP. Advances in physiological and molecular aspects of plant cold tolerance. J Plant Interact. 2017;12(1):143-57. https://doi.org/10.1080/17429145.2017.1308568.

10. Stockinger E, Gilmour S, Thomashow M. *Arabidopsis thaliana* CBF1 encodes an AP2 domain-containing transcriptional activator that binds to the C-repeat/DRE, a cis-acting DNA regulatory element that stimulates
**Additional file 13: Table S7.** Description of accessions of the pea reference collection for their haplotype at the frost damage (FD)-associated loci of the GWA study and at the *Hr* locus. At the linkage groups (LGs) I, II, III, V, VI and VII, the favorable haplotypes are shown by a red background, the unfavorable haplotypes by a blue background, as described in the legend of Additional file 10. Accessions with undefined haplotypes or intermediate haplotypes were not presented. The same colour code has been used to describe the favorable (red background: *Hr*) and unfavorable (blue background: *hr*) allele for the *Hr* gene, as determined in Lejeune-Hénaut et al. [21]. The mean frost damage score observed in the field experiment as well as its standard error (SE) are given. Frost scores are ranging from 0 (no damage) to 5 (dead plant). The passport data of accessions are extracted from the Additional file 6.

**Additional file 14: Table S8.** Daily air and soil temperatures during the field experiment in Clermont-Ferrand Theix.

## Abbreviations

AP2/EREBP: APETALA2/Ethylene Responsive Element Binding Protein;
BIC: Bayesian Information Criterion; BLUPS: Best Linear Unbiased Predictions;
CBF: C-repeat Binding Factor; cM: centi Morgan; COR: Cold Regulated;
CVg: Coefficient of Variation of the genetic variance; DAPC: Discriminant
Analysis of Principal Components; DREB: Dehydration-Responsive Element-
Binding; Elf3: Early flowering 3; FD: Frost Damage; FD_CC: Frost damages in
the controlled conditions experiment; FD_Field_date1: Frost damages in the
field experiment at the first date of scoring; FD_Field_date2: Frost damages
in the field experiment at the second date of scoring; FD_Field_date3: Frost
damages in the field experiment at the third date of scoring;
FD_Field_date4: Frost damages in the field experiment at the fourth date of
scoring; GWAS: Genome Wide Association Study; H²: Broad sense heritability;
Hr: High response to photoperiod; INRAE: National Research Institute for
Agriculture, Food and Environment; k: DAPC Cluster; K: Kinship matrix;
K_LG: Kinship matrix corresponding to a given Linkage Group; LD: Linkage
Disequilibrium; LG: Linkage Group; LMM: Linear Mixed Model; MAF: Minor
Allele Frequency; NAM: Nested Association Mapping; NIL: Near Isogenic Line;
PCs: Principal Components; PCA: Principal Components Analysis; Q: Structure
matrix; QTL: Quantitative Trait Locus; RIL: Recombinant Inbred Line;
SE: Standard Error; SNK: Student-Newman-Keuls; SNP: Single Nucleotide
Polymorphism; SSD: Single Seed Descent; SSR: Simple Sequence Repeat;
Vg: Genetic variance; WFD: Winter Frost Damage

## Acknowledgements

The authors would like to thank the INRAE experimental units of Estrées-Mons, Dijon and Theirx, France, for their contribution to the field experiments. They are also grateful to Catherine Desmetz (Agroécologie, INRAE) for additional genotyping of the reference collection, to Marianne Chabert-Martinello, Anthony Klein and Jean-Bernard Magnin-Robert (Agroécologie, INRAE) for providing informations about the constitution of the reference collection and to Renaud Rincent (GDEC, INRAE) for his advices in genome wide association analyses.

## Authors' contributions

VF, RD and MT collected phenotypic data in the field and controlled conditions experiments. GA and JB coordinated the genotyping of the association mapping collection and provided genotyping data. SB performed the phenotypic and genetic analyses, as well as association mapping and haplotypes analysis and wrote the manuscript. MS and SSN assisted in analysing genetic data by testing GWAS models and providing a part of R scripts, respectively. ILH and BD conceived and coordinated the study. ILH managed the technical and financial supplies of the study within the PeaMUST project. JLH, BD and ILH coordinated the overall progress and financial support of the study. ILH, BD and NB reviewed and contributed to draft the manuscript. All authors read and approved the final manuscript.

## Funding

This work was supported by the PhD fellowship of Sana Beji which was co-funded by the University of Lille, doctoral school for materials, radiation and environmental sciences (France) and by the region Hauts-de-France (France). The acquisition of the phenotyping data was funded by the SNPEA project (ANR06-GPLA019). The acquisition of the genotyping data was funded by the GENOPEA project (ANR-09-GENM-026). The analyzes, interpretation of

data and writing of the manuscript were realized in the frame of the Pea-MUST project (ANR11-BTBR-0002). The three project fundings come from the French government and were managed by the Research National Agency (ANR). The funding agency did not contribute to the design of the study or collection, nor to the analysis and interpretation of data or to the writing of the manuscript.

## Availability of data and materials

Genetic positions and sequences of markers used for the current study are available in Tayeh et al. [36]. The pea genome assembly v.1a explored in this study is available at https://urgi.versailles.inra.fr/Species/Pisum (Kreplak et al. [39]). The data supporting the findings of this study were used under license within the PeaMUST project and are not publicly available. Phenotyping and genotyping data, all other intermediate data files and scripts are however available from the authors upon reasonable request and subjected to data transfer agreement.

## Ethics approval and consent to participate

Not applicable.

## Consent for publication

Not applicable.

## Competing interests

The authors declare that they have no competing interests.

## Author details

¹BioEcoAgro, INRAE, Univ. Liège, Univ. Lille, Univ. Picardie Jules Verne, 2, Chaussée Brunehaut, F-80203 Estrées-Mons, France. ²GCIE-Piccardie, INRAE, F-80203 Estrées-Mons, France. ³GQE - Le Moulon, INRAE, Univ. Paris-Sud, CNRS, AgroParisTech, Univ. Paris-Saclay, F-91190 Gif-sur-Yvette, France.

⁴Agroécologie, AgroSup Dijon, INRAE, Univ. Bourgogne, Univ. Bourgogne Franche-Comté, F-21000 Dijon, France.

Received: 9 August 2019 Accepted: 20 July 2020

Published online: 04 August 2020

## References

1. FAOSTAT. FAOSTAT data base. Food and agriculture organization of the United Nations - Statistic Division. 2018. http://www.fao.org/faostat/en/#data/. Accessed 27 Mar 2020.

2. Graham P, Vance C. Legumes: importance and constraints to greater use. Plant Physiol. 2003;131:872-7. https://doi.org/10.1104/pp.017004.

3. Jeuffroy MH, Baranger E, Carrouee B, de Chezelles E, Gosme M, Henault C, et al. Nitrous oxide emissions from crop rotations including wheat, oilseed rape and dry peas. Biogeosciences. 2013;10(3):1787-97. https://doi.org/10.5194/bg-10-1787-2013.

4. Benezit M, Biarnes V, Jeuffroy MH. Impact of climate and diseases on pea yields: what perspectives with climate change? OCL. 2017;24(1):D103. https://doi.org/10.1051/ocl/2016055.

5. Lejeune-Hénaut I, Bourion V, Eteve G, Cunot E, Delhaye K, Desmyter C. Floral initiation in field-grown forage peas is delayed to a greater extent by short photoperiods, than in other types of European varieties. Euphytica. 1999; 109(3):201-11. https://doi.org/10.1023/a:1003727324475.

6. Fowler DB, Breton G, Limin AE, Mahfoozi S, Sarhan F. Photoperiod and temperature interactions regulate low-temperature-induced gene expression in barley. Plant Physiol. 2001;127(4):1676-81. https://doi.org/10.1104/pp.010483.

7. Thomashow M. Plant cold acclimation, freezing tolerance genes and regulatory mechanisms. Annu Rev Plant Physiol Plant Mol Biol. 1999;50:571-99. https://doi.org/10.1146/annurev.arplant.50.1.571.

8. Ruelland E, Vaultier M, Zachowski A, Hurry V. Cold signalling and cold acclimation in plants. Adv Bot Res. 2009;49:35-146. https://doi.org/10.1016/S0065-2296(08)00602-2.

9. Rihan HZ, Al-Issawi M, Fuller MP. Advances in physiological and molecular aspects of plant cold tolerance. J Plant Interact. 2017;12(1):143-57. https://doi.org/10.1080/17429145.2017.1308568.

10. Stockinger E, Gilmour S, Thomashow M. *Arabidopsis thaliana* CBF1 encodes an AP2 domain-containing transcriptional activator that binds to the C-repeat/DRE, a cis-acting DNA regulatory element that stimulates