Genetic Dissection of Heat Stress Tolerance in Faba Bean (Vicia faba L.) Using GWAS

Heat waves are expected to become more frequent and intense, which will impact faba bean cultivation globally. Conventional breeding methods are effective but take considerable time to achieve breeding goals, and, therefore, the identification of molecular markers associated with key genes controlling heat tolerance can facilitate and accelerate efficient variety development. We phenotyped 134 accessions in six open field experiments during summer seasons at Terbol, Lebanon, at Hudeiba, Sudan, and at Central Ferry, WA, USA from 2015 to 2018. These accessions were genotyped using genotyping by sequencing (GBS), and 10,794 high quality single nucleotide polymorphisms (SNPs) were discovered. These accessions were clustered in one diverse large group, although several discrete groups may exist surrounding it. Fifteen lines belonging to different botanical groups were identified as tolerant to heat. SNPs associated with heat tolerance using single-trait (ST) and multi-trait (MT) genome-wide association studies (GWASs) showed 9 and 11 significant associations, respectively. Through the annotation of the discovered significant SNPs, we found that SNPs from transcription factor helix–loop–helix bHLH143-like S-adenosylmethionine carrier, putative pentatricopeptide repeat-containing protein At5g08310, protein NLP8-like, and photosystem II reaction center PSB28 proteins are associated with heat tolerance.


Introduction
Faba bean (Vicia faba L.) is the fourth most important cool season legume crop grown globally in diverse cropping systems and environments [1]. It plays a valuable role in increasing food production and enhancing sustainable farming methods due to its ability to fix nitrogen and improve soil structure [2,3]. Faba bean is currently cultivated on 2.65 million ha worldwide, with a total dry grains production of 5.6 Mt [4]. This crop is used as human food and animal feed, mainly for pigs, horses, poultry, and pigeons, in many countries [5] because it is rich in protein and essential amino acids [6] such as arginine, lysine, and leucine. It also has high iron (22 to 78 mg/kg) and zinc (45-61 mg/kg) contents [7]. On the other hand, several constraints, such as biotic and abiotic stresses, significantly affect its productivity. Heat stress is one of the major constraints affecting faba bean in heat prone environments.
Due to climate change, heat waves (periods of abnormally hot weather lasting from days to weeks) are expected to become more intense and frequent, and the average global temperature is also expected to increase by 1.5 • C, as reported by the assessments of the Intergovernmental Panel on Climate Change (IPCC) (https://www.ipcc.ch/sr15/, accessed on 5 July 2021). This will have a severe impact on crop production in many regions, especially in the Mediterranean region and Nile Valley countries such as Ethiopia, Egypt, and Sudan, where faba bean is grown extensively as a major food legume crop. Moreover, the projected climate change will affect the stability of faba bean cultivars, which highlights the need to develop heat tolerant faba bean cultivars with stable floral and reproductive development [8]. Heat stress can significantly reduce faba bean yield [9], especially when the temperature rises above 30 • C [10]. Heat stress during flowering reduces pollen viability, which reduces the pod set [11]. The effects of heat stress on pollen development have also been reported in many other legumes such as chickpea (Cicer arietinum) [11], common bean (Phaseolus vulgaris) [12], cowpea (Vigna unguiculata) [13], groundnut (Arachis hypogaea) [14], soybean (Glycime max) [15], and lentil (Lens culinaris) [16].
The field selection method is used in heat prone environments to identify sources of heat tolerance and results in the development of heat tolerant varieties [7]. These conventional breeding tools are effective but time consuming. Although recent advances in the next generation sequencing technologies (NGS) have facilitated the generation of large volumes of sequences [17][18][19], the large genome size of faba bean (~13 Gb) has so far limited the development of comprehensive genomic information [20][21][22]. However, NGS can be easily used to discover single nucleotide polymorphisms (SNPs) associated with key traits through genome-wide association studies (GWASs) [17,22,23]. GWASs utilize the concept of linkage disequilibrium (LD) in diverse populations to detect SNPs associated with different traits [24]. One of the key factors that affects the power of GWASs is the size and structure of the population used for the analysis [25]. The population can be computationally increased using multi-variate (environments/traits) genome-wide association study (M-GWAS) models if the evaluation trial is not large enough [26].
Breeding new genotypes with tolerance/avoidance to projected stresses is important when it comes to thoroughly quantifying the crop's response to heat stress [27]. So far, efforts made to explore the physiological and molecular basis of heat tolerance have been mainly focused on vegetative growth [28] and reproductive stages [11,29]. SNPs associated with traits of interest in faba bean have been reported for some biotic stresses such as Ascochyta blight resistance [30,31] and broomrape resistance [31], for quality traits including low vicine-convicine content [32], and for herbicide tolerance [33]. However, no comprehensive study has been conducted to detect molecular markers associated with heat tolerance.
The major objective of the present research was to identify markers associated with heat tolerance in a diverse faba bean collection. The secondary objective was to study the relationship between grain yield and different physiological and agronomical traits in diverse heat prone environments.

Phenological and Yield Traits
Combined variance analysis using residual maximum likelihood (REML) and spatial models for each trial are presented in Table 1. Significant differences were observed among genotypes and among environments for all evaluated traits. The genotype × environment (G × E) interaction (Table 2) were significant for all traits except for grain yield, indicating that, contrary to the order of other traits, the order of yield performance of the tested accessions remained the same across the environments. Spatial analysis for each season-location combination indicated a wide variability among genotypes for days to flowering (DFLR) and days to maturity (DMAT) at Terbol and Hudeiba, from 2015 to 2017. Significant differences were observed among genotypes for grain yield (p < 0.001 at all locations/seasons in 2016 and 2017 and p < 0.05 at Terbol in 2015). Significant differences were also found among genotypes for the number of seeds per plant and hundred seed weight in all environments. Plant height had significant variation among genotypes at all locations except for Hudeiba in 2017. Means, ranges, and standard errors scored for all traits and locations are presented in Table 3. The average grain yield varied from 1.3 g at Hudeiba station to 13.9 g at the experimental station in Central Ferry, indicating wide variability across locations. Large variability between locations were also reported for other traits. The hundred seed weight showed consistently lower values under heat stress at Hudeiba, Sudan.

Physiological Traits
Combined variance analysis (Table 1) showed significant differences among accessions for PG and CL. However, no significant difference was observed for canopy temperature (CT) in both seasons at Hudeiba station, Sudan. Significant differences were detected among environments (p < 0.001), and accession X environment were significant for PG and CL. Means and ranges for CL and PG are presented in Table 3. The average PG varied from 8.87% in 2015 to 28.4% at Terbol station, while CL varied from 40.3 to 41.2% at Hudeiba.

Relationship between Traits
Since CT did not vary significantly among accessions, it was excluded from the correlation and regressions analysis. Correlations and stepwise regression analysis using mean estimates are presented in Table 4. GY was significantly and positively correlated with NPP (p < 0.001), PG and PLHT (p < 0.01), and NSP (p < 0.05). HSW was significantly and positively correlated with PLHT (p < 0.001) but negatively correlated with PG (p < 0.001). CL was highly correlated with PLHT, but it was not correlated with GYP under heat stress. Therefore, CL appeared inefficient for selection for yield in a heat prone environment, despite it being an effective indicator of plant biomass. Table 4. Correlation analysis between mean phenotyping estimates for different phenotyping and physiological traits and regression equation between grain yield and associated traits.  DFLR was positively correlated with DMAT and negatively correlated with HSW and PLHT. Early flowering lines with late maturity showed higher seeds sizes than late flowering lines. The regression coefficient (Table 4) indicated that GYP had moderate correlation with NPP (α = 0.31), PG (α = 0.1), and HSW (α = 0.01). These results indicated that PG and NPP are effective parameters to select for heat tolerance in faba bean. A total of 15 lines from different subspecies ( Table 5) were selected based on grain yield, pollen germination, and plant height, of which eight were equina types, two belong to major types, three to minor types and two to paucijuga.

Population Structure
A principal component analysis (PCA) is presented in Figure 1. The results showed that the tested accessions were clustered in one large group, although three discrete groups of genotypes may exist surrounding it. Labelling each faba bean genotype according to its subspecies or geographic origin does not discriminate any clear group, suggesting that the assembled faba bean germplasm is composed of a single divergent population ( Figure 1).

Genome-Wide Association Analysis
Analysis of 39 traits across six environments using ST-GWAS revealed 68 SNP-trait associations for 25 traits, while the remaining 14 traits had no associations. Figure 2 depicts the QQ plots that were used to compare observed vs. expected values to demonstrate how well the GWAS model performed in identifying trait-associated SNPs. All traits except CL exhibited associations in at least one environment. The highest number of associations were detected for NSP, NPP, and GYP in Terbol summer 2016 with 9, 7, and 6 associations, respectively. Out of these 68 associations, only 14 were significant following Bonferroni correction at p < 4.6 × 10 −6 while the remaining 54 associations were suggestive (Supplementary Table 1). Nine SNPs showed associations with multiple traits/environments, of which the SNP SCONTIG60075_82 showed the highest number of four associations with GYP, NPP, and NSP in Hudeiba 2017 as well as NSP in Hubeiba 2016. The SNP SNODE_39949_LENGTH_82_COV_1.024390_37 showed associations with NPP, GYP, and NSP in Hudeiba 2017. The remaining eight SNPs showed associations with two traits, with five being associated with NSP and NPP in Terbol summer 2016, two associated with GYP and NSP in Hudeiba 2016, and the remaining one associated with NPP and GYP in Terbol summer 2015 (Supplementary Table S1). The MT-GWAS model which fit each of the 10 studied traits across all environments in a single analysis detected a total of 56 SNP-trait associations, of which 11 were significant when considering Bonferroni correction (Table 6, Figure 2). These associations were clustered in 50 QTL based on linkage disequilibrium (LD) analysis. NPP showed 10 QTL, followed by NSP, which had nine QTL, while DMAT and CL showed one QTL each. Each of the following traits, GYP, HSW, and PG, had seven QTL, while PLHT showed four QTL and CT and DFLR had only two QTL each (Supplementary Table S2). Like the ST-GWAS analysis, the SNPs SNODE_39949_LENGTH_82_COV_1.024390_37 and SCONTIG60075_82 showed association with three traits, GYP, NPP, and NSP. The SNP SNODE_40333_LENGTH_77_COV_34.987015_87 was associated with GYP and NPP, while the SNP SCONTIG124089_41 was associated with both NPP and NSP.

Genome-Wide Association Analysis
Analysis of 39 traits across six environments using ST-GWAS revealed 68 SNP-trait associations for 25 traits, while the remaining 14 traits had no associations. Figure 2 depicts the QQ plots that were used to compare observed vs. expected values to demonstrate how well the GWAS model performed in identifying trait-associated SNPs. All traits except CL exhibited associations in at least one environment. The highest number of associations were detected for NSP, NPP, and GYP in Terbol summer 2016 with 9, 7, and 6 associations, respectively. Out of these 68 associations, only 14 were significant following Bonferroni correction at p < 4.6 × 10 −6 while the remaining 54 associations were suggestive (Supplementary Table S1). Nine SNPs showed associations with multiple traits/environments, of which the SNP SCONTIG60075_82 showed the highest number of four associations with GYP, NPP, and NSP in Hudeiba 2017 as well as NSP in Hubeiba 2016. The SNP SNODE_39949_LENGTH_82_COV_1.024390_37 showed associations with NPP, GYP, and NSP in Hudeiba 2017. The remaining eight SNPs showed associations with two traits, with five being associated with NSP and NPP in Terbol summer 2016, two associated with GYP and NSP in Hudeiba 2016, and the remaining one associated with NPP and GYP in Terbol summer 2015 (Supplementary Table S1).   The MT-GWAS model which fit each of the 10 studied traits across all environments in a single analysis detected a total of 56 SNP-trait associations, of which 11 were significant when considering Bonferroni correction (Table 6, Figure 2). These associations were clustered in 50 QTL based on linkage disequilibrium (LD) analysis. NPP showed 10 QTL, followed by NSP, which had nine QTL, while DMAT and CL showed one QTL each. Each of the following traits, GYP, HSW, and PG, had seven QTL, while PLHT showed four QTL and CT and DFLR had only two QTL each (Supplementary Table S2). Like the ST-GWAS analysis, the SNPs SNODE_39949_LENGTH_82_COV_1.024390_37 and SCONTIG60075_82 showed association with three traits, GYP, NPP, and NSP. The SNP SNODE_40333_LENGTH_77_COV_34.987015_87 was associated with GYP and NPP, while the SNP SCONTIG124089_41 was associated with both NPP and NSP.

Candidate Gene Annotation
Ten identified SNPs located within diverse functional genes are presented in Tables 7 and 8. These SNP markers were associated with PG, HSW, NSP, and GYP using MT-GWAS. The annotated SNP markers associated with PG were: (1) Contig82855 located within one of two probable polygalacturonase genes related to lupin (Lupinus angustifolius), jequirity bean (Abrus precatorius), and the pseudomolecule Pd05 gene of sweet almond (Prunus dulcis); (2) NODE_38942_length_69_cov_118.768112 located within one of seven uncharacterized LOC100810394 genes related to soybean, jequirity bean, chickpea, and barrel clover (Medicago truncatula); (3) NODE_6662_length_69_cov_474.000000 located within one of three genes of transcription factor-basic helix-loop-helix (TF-bHLH143) related to chickpea, barrel clover, and lupin; (4) NODE_7398_length_62_cov_214.516129 located within six genes of S-adenosylmethionine carrier 1 related to chickpea and barrel clover; and (5) NODE_7979_length_116_cov_512.344849 located within one of seven putative pentatricopeptide repeat-containing protein At5g08310 genes related to chickpea and barrel clover. Table 7. Gene annotations of different identified SNPs markers associated with pollen germination.

Markers
Gene Names
In addition, SNP Contig82391 and NODE_14795_length_67_cov_68.791046 were respectively associated with NSP located within one of two proteins, NLP8-like (LOC101496898) genes of chickpea, and with GYP located within the LOC25500962 gene of barrel clover.
Finally, Contig60075, which was a major SNP marker identified by ST-GWAS, was located within photosystem II reaction center PSB28 protein. This marker is associated with DFLR, DMAT SNP, NPP, and GYP.

Discussion
Heat stress during the reproductive phase significantly affects faba bean yield. The identification of a heat tolerant faba bean is therefore required to assure appropriate yield in heat prone environments. Yield and yield component-related traits are highly affected by heat stress in faba bean [8,10]. The results of our field experiments confirmed the results of earlier studies that heat stress negatively affects faba bean yield and yield components such as pollen germination, pod number, seed number, and grain yield. Grain yield was found to be positively associated with pod and seed numbers. Heat stress therefore damages faba bean during the flowering period, when the viability of pollen is critical for successful reproduction [8]. Moreover, our experiments were conducted under heat prone and wellwatered conditions, and therefore canopy temperature was not associated with grain yield, as reported in wheat [34]. Moreover, in agreement with the results obtained under water stress conditions in faba bean [35], chlorophyll content and plant height were not associated with grain yield.
To date, conventional breeding methods have identified only limited heat tolerant faba bean cultivars [36]. These methods are relatively slow in achieving significant progress. The development of genetic markers associated with important traits under heat stress in faba bean is essential to improve selection accuracy and make selections in early generations. Such markers are effective for selecting adapted and tolerant accessions. Few published resources have used SNP genotyping to provide a useful genetic map that could be applied in faba bean breeding [17]. In the present study, GBS technology was used to genotype 134 faba bean accessions that belong to four different botanical varieties originated from either the collection of landraces from 42 countries or from ICARDA breeding program. Our GBS pipeline produced 10,794 SNPs on a diverse faba bean germplasm. Recently, SNP genotyping analysis was used to generate a collection of 1819 gene-based SNP markers in three recombinant inbred line (RIL) populations to enable the design of a high-density consensus genetic map in faba bean [22]. The high number of SNPs produced in the present study may be due to the high heterogeneity, as we used diverse genotypes. A parental population with a common ancestor produces a low amount of genetic variability and therefore a low number of SNPs [37].
Despite the large genetic variability observed in the faba bean germplasm, only one large heterogeneous group was observed by PCA (Figure 1), confirming the lack of genetic or sterility barriers between different subspecies [38]. This contrasts with the previous ADMIXTURE analysis that reported the presence of two ancestral subpopulations within the germplasm [33]. While the germplasm used in both studies largely overlap, the previous one included 52 additional accessions, all belonging to the first subpopulation, while the shared accessions all belonged to the second subpopulation.
A limited number of studies have investigated the genetic control of different traits in faba bean through biparental populations, multi-parental populations, or GWAS [22,30,33,39], but none of these targeted QTL associated with heat tolerance. Thus, the QTL identified for heat tolerance in the present study would be very useful in faba bean breeding programs. Moreover, the QTL detected with the MT-GWAS model are supposed to be stable across environments, giving them a broader practical impact in breeding [26,40] Heat stress has harmful effects on plant metabolism and impacts crop growth and development throughout vegetative and reproductive phases [41]. Plants exposed to heat stress generate reactive oxygen species (ROS), antioxidant substances, plant hormones, and other secondary metabolites [41]. Therefore, heat stress response pathways rely heavily on transcription factor-mediated gene expression regulatory networks. Annotated genes related to SNPs under heat stress are the factors regulating heat tolerance in faba bean. As an example, SNPs associated with PG under heat stress were located within transcription factor (TF) bHLH143-like, probable polygalacturonate inhibiting protein (PGIP), S-adenosylmethionine carrier 1 (SAMC), and putative pentatricopeptide repeat-containing protein At5g08310 (RCP). Therefore, these genes may have crucial roles in controlling pollen germination (PG) under heat stress. Similar results showed that TF and PPRCP genes had crucial roles in the resistance of Benincasa hispida to high temperature [42]. In addition, SMAC was considered as an important precursor of polyamine and ethylene and leads to stress tolerance in plants [43].
Psb28 is an extrinsic protein of photosystem II (PSII) that is found in photosynthetic organisms ranging from cyanobacteria to higher plants. Several studies have revealed that Psb28 is involved in the recovery of PSII under high temperature stress [44]. Additionally, the correlation between gene expression of Psb28 and plant response to heat stress was detected in several plant species including Populus tomentosa and perennial ryegrass [45,46] The SNPs marker Contig82391 located with Cicer arietinum protein NLP8-like (LOC101496898) was found to be associated with the number of seeds per plant under heat stress. NLP8 was described in Arabidopsis as an essential protein for nitrate-promoted seed germination [47].

Genetic Materials
The present study involved a set of 134 faba bean accessions comprising 16 breeding lines developed at ICARDA and 118 landraces originating from 42 countries. The sub-species of the germplasm set was characterized as 42 major, 17 minor, 66 equina, and 9 paucijuga. This germplasm set represents wide diversity, as it represents the 1000 accessions previously assessed by Maalouf et al. [20].

Experimental Designs
Six experiments were conducted in an incomplete block design with two replicates at three locations between 2015 and 2017 (Table 9). Experiments details are described as follows: (1) Three experiments were conducted during the summer season at the ICARDA Terbol station (35.98 • N, 33.88 • E, 890 m above sea level). The summer season at the Terbol station, Lebanon runs from June to October and is characterized by hot and dry weather, with temperatures above 35 • C during the flowering and pod set time (Figure 3a). A 50 mm irrigation was provided to the crop every seven days using drip irrigations to ensure enough moisture for crop growth. The soil in Terbol station is a deep and rich clay loam soil.

Phenotyping for Heat Stress
The population was phenotyped for days to 50% flowering (DFLR) and days to 90% maturity (DMAT) on plot basis and plant height (PLHT), number of pods per plant (NPP), number of seeds per plant (NSP), and grain yield per plant (GYP) as an average of three  The plot size in all six experiments was 2 m in length with two rows 45 cm apart. Soil fertilization of NPK 15-15-15 at 250 kg ha −1 was added; weeds were controlled by the pre-emergence application of pendimethalin at 1200 g a.i./ha followed by manual weeding; insects and fungal diseases were controlled regularly to eliminate any additional stress on plants other than the heat stress.

Phenotyping for Heat Stress
The population was phenotyped for days to 50% flowering (DFLR) and days to 90% maturity (DMAT) on plot basis and plant height (PLHT), number of pods per plant (NPP), number of seeds per plant (NSP), and grain yield per plant (GYP) as an average of three random plants per plot. Hundred seed weight (HSW) was recorded as the weight of a hundred random seeds per plot. In addition, observations on chlorophyll content and canopy temperature were also recorded at Hudeiba research station during both seasons.
Chlorophyll content (Cl): The chlorophyll content of three random leaves per plot was measured during morning hours at the end of the flowering period using a chlorophyll meter (SPAD-502 Minolta (Spectrum Technologies Inc., Plainfield, IL, USA). The SPAD 502 Plus Chlorophyll Meter instantly measures chlorophyll contents to assess the effect of heat stress on vegetative growth and its relation to yield and other agronomic traits.
Pollen Germination (PG): The pollen germination test was conducted during three consecutive seasons from 2015 to 2017 at Terbol station. Fresh pollens were randomly collected from 15 opened flowers of each plot when the day temperature exceeded 32 • C. The collected pollens were then kept in microfuge tubes in iceboxes for a few hours prior to the germination test. The pollen germination test was conducted in two replicates per plot following the method described earlier [48]. A medium containing 0.1 g/L boric acid (H3BO3), 0.3 g/L calcium nitrate (Ca (NO3)2), 100 g/L sucrose, and 10 g/L agar (Sigma, Aldrich-Germany) was prepared and autoclaved at 121 • C before being poured into Petri dishes (20 mL/dish of 9 cm diameter). The pH was adjusted to 6.3. The freshly collected pollen grains were dusted and dispersed on the germination surface using a needle and then incubated at 37 • C for 16 h before the germination was terminated by adding a drop of acetocarmine to the medium. Pollen grains were counted as germinated if the elongation of pollen tubes were at least equal to the diameter of the grain. The germination rate was determined by counting the germinated pollens per hundred observed pollen grains in each replicate under a Cole-Parmer DX53056101 Stereozoom Microscope, with a built-in 2.0 mega pixel digital camera, 110-220 VAC.

DNA Extraction and Genome by Sequencing Analysis
Details relating to the DNA extraction, GBS protocol, and data analysis were previously described in Abou Khater et al. [33]. In brief, genomic DNA was extracted using the Qiagen Plant DNA Preparation Kit, while the two restriction enzymes used to generate the GBS libraries, PstI and MspI, were prepared with 48 barcode adapters with a 4-9 bp sequence [49] and produced approximately four million single reads (100 bp) per genotype using the Illumina HiSeq 2500 (San Diego, CA, USA) using standard Illumina protocols and kits which produce high-output paired-end 100 bp reads (Montreal, QC, Canada). The TASSEL-GBS 3.0 bioinformatics analysis package was used to process raw sequence files in FASTQ format using the first steps of the TASSEL-GBS pipeline with its default parameters. The TASSEL-GBS workflow was used to trim reads, remove barcodes, and generate a single tag count file for each sample. The raw sequence files in FASTQ format were assembled and converted into FASTA format, from which a faba bean reference genome was constructed using NCBI and pulsedb genomic and transcriptomic sequences, as well as our newly assembled sequences. This step should result in a unified reference genome that covers more transcriptionally active regions of the faba genome. GBS tags were aligned to the generated reference genome using Bowtie2 V2.2.4 [50] with the very-sensitive-local option. The generated BAM files were used in the TASSEL-GBS pipeline for SNP marker calling. SNPs with >85% call rate and >5% minor allele frequency were used in subsequent analyses.

Statistical Analysis of Phenotyping Data
Variance analysis using the spatial model for all traits was applied using GenStat V19 with the automatic incomplete block design modules [51]. The replicates and block within replicates were fitted as random variables and accessions were considered as fixed variables. In addition to variance parameters, the output of the above model displays Wald statistics, the predicted mean value of each accession, and the standard error. Differences among accessions and environments were examined using the multi-environment trials analysis (META), in which accessions were fitted as fixed variables while environments and the accession × environment interaction (G × E) were fitted as random variables. The META analysis was conducted using the method of residual maximum likelihood (REML).
Correlation analysis and step-wise regression analysis were conducted to examine the relationship between traits and to define traits that have an influence on grain yield in faba bean under heat prone conditions. The combined mean estimates were obtained and used for correlations and step-wise multiple regression analyses. The selected heat tolerant accessions are supposed to have higher yield and higher values for traits influencing yield under stress conditions.

Population Structure and Genome-Wide Association Analysis
PCA was estimated using TASSEL 5 [52] on the 10,794 SNPs. The PCA vectors were then plotted in R [53] with the ggplot2 R package [54]. Faba bean accessions were labelled with different colors and shapes according to their geographic origin and subspecies, respectively.
Two GWAS models were used in the present study, the single-trait (ST) model, in which each environment is fitted independently, and the multi-trait (MT) model, in which all environments are jointly fitted in a single analysis. GEMMA software was used to fit both models with its default parameters and was also used to calculate the genomic relatedness matrix (GRM) following the VanRaden (2008) method [55]. Significant threshold was defined with Bonferroni correction at p < 0.05 and all SNP-trait associations that had p < 1 × 10 −4 were reported as suggestive associations. To determine if different significant SNPs for each analysis were associated with the same quantitative trait locus (QTL) or different QTL, the linkage disequilibrium (LD) among them was calculated following the Weir method [56]. The sequences flanking associated SNPs were blasted against the NCBI database to annotate potential candidate genes underlying the causal variants.

Conclusions
Heat stress is one of the major constraints affecting faba bean in many production regions, especially when waves of high temperatures coincide with the reproductive phase. It also limits the crops expansion in new geographic regions where high temperatures are common during the growing season. By improving heat tolerance in faba bean, the crop can be expanded to new niches in sub-Saharan Africa and therefore contribute to the improvement of soil quality. The identified heat tolerant lines can be used to develop cultivars adapted to heat prone areas and potentially contribute to enhanced food security in countries such as Sudan, Egypt, and Ethiopia. Considering the importance of heat tolerance in faba bean, it is imperative to breed elite cultivars that feature this trait. However, field selection is very costly and time consuming, and this is reflected by low genetic gains over the years. The integration of genomic selection and marker-assisted selection will improve selection accuracy and intensity and will shorten the breeding cycle when selecting at early generations. In the present study, we identified SNP markers associated with agronomic and physiological traits under heat prone conditions using single-and multi-trait association. These SNP markers located within functional genes expressed under heat stress will facilitate and fasten the development of heat tolerant cultivars and can be used to achieve the introgression of heat tolerant genes into a desired agronomic background.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11091108/s1, Table S1: SNP associated with different traits under heat stress at Terbol summer (TRS) and Hudeiba (HD) in 2016 and 2017 seasons using ST-GWAS model; Table S2: Multi-trait genome-wide association study (MT-GWAS) results for grain yield and yield components. Underlined p-values are the significant p-values while underlined SNPs are the SNPs that were detected in multiple traits.  Data Availability Statement: The datasets generated and analyzed during the current study are available from the corresponding authors on request.