Genome-Wide Association Mapping for Heat Stress Responsive Traits in Field Pea

Environmental stress hampers pea productivity. To understand the genetic basis of heat resistance, a genome-wide association study (GWAS) was conducted on six stress responsive traits of physiological and agronomic importance in pea, with an objective to identify the genetic loci associated with these traits. One hundred and thirty-five genetically diverse pea accessions from major pea growing areas of the world were phenotyped in field trials across five environments, under generally ambient (control) and heat stress conditions. Statistical analysis of phenotype indicated significant effects of genotype (G), environment (E), and G × E interaction for all traits. A total of 16,877 known high-quality SNPs were used for association analysis to determine marker-trait associations (MTA). We identified 32 MTAs that were consistent in at least three environments for association with the traits of stress resistance: six for chlorophyll concentration measured by a soil plant analysis development meter; two each for photochemical reflectance index and canopy temperature; seven for reproductive stem length; six for internode length; and nine for pod number. Forty-eight candidate genes were identified within 15 kb distance of these markers. The identified markers and candidate genes have potential for marker-assisted selection towards the development of heat resistant pea cultivars.


Introduction
Pea (Pisum sativum L., 2n = 14) is a major pulse crop widely grown in the temperate regions primarily for its nutritional values as a source of protein, slowly digestible starch, essential minerals, high fiber and low fat; and soil fertility benefits as it fixes atmospheric nitrogen [1][2][3]. However, as a cool season crop, pea is prone to heat and drought stress, with warm summers causing shortened life cycles, abortion of floral components and pods, and thus economic yield loss [4][5][6]. Due to global warming, the average surface temperature is predicted to increase by 3.7 • C by the end of this century, and thus heat stress is expected to be even more challenging in the future [7].
Genetic improvement of pea for heat and drought resistance is a promising approach to stabilize yield under environmental stresses. Pea germplasm has a wide range of diversity in morpho-anatomical, biochemical and physiological characteristics [8,9]. Among other things, such diversity has been explored to identify traits associated with heat response [10][11][12]. Pigments including chlorophylls, carotenoids, anthocyanins contribute to heat tolerance through heat dissipation and protection of season [5]. Impact of heat and drought is severe when it occurs during reproductive stages. Saskatoon 2015 was the most stressed environment as indicated by mean daily maximum air temperatures > 27 • C, 18 days where air temperature was > 28 • C, and drier conditions during the reproductive stage. Similarly, 2017 Saskatoon was also under heat and drought stress during the reproductive stage with average air temperature~26 • C, 16 days where air temperature was > 28 • C, and relatively low total precipitation. The remaining three environments were generally ambient and considered as control environments (Table 1). Table 1. Seeding date, average maximum, minimum and 24 h daily mean temperatures, number of days when the daily maximum temperature was greater than 28 • C, and total monthly precipitation at different growth and development stages of pea at each environment. Note: temp, temperature; mm, millimeter; Means of environmental variables that do not share a letter within a column under each growth stage are significantly different from each other. The data were analyzed by one-way ANOVA and followed with the Tukey-HSD test for the mean separations.

Phenotypic Measurements, Analysis of Variance, and Marker Detection through Association Mapping
Variance components of genotype (G), environment (E), and G × E interaction together with their significance on the six traits used in this study is presented in Table 2. For all traits analyzed, normality of residuals and homogeneity of variance were met. Note: * Significant at the 0.05 level of probability; ** Significant at the 0.01 level of probability; *** Significant at the 0.001 level of probability; ns, not significant at the 0.05 level. SPAD, soil plant analysis development; PRI, photochemical reflectance index.
Descriptive statistics for minimum, maximum and mean values of phenotypic measurements on the traits of the GWAS panel across five environments are summarized in Table 3 and Figure 1. Chlorophyll concentration, measured by a SPAD meter, was affected by genotype, environment and their interaction; and the variance component analysis showed that maximum variation (67.9%) among the GWAS panel was due to the genotype effect, and the broad sense heritability was 0.95. Overall, genotype chlorophyll concentration ranged from 26.6 to 57.6 SPAD values under heat stress, and 30.0 to 67.5 under control conditions (Table 3). On average, the heat stressed environments had 3% less SPAD value than the ambient environments. Six markers (Chr5LG3_150942510, Chr5LG3_446272814, Chr5LG3_449362407, Chr5LG3_566189589, Chr5LG3_569788697, and Chr5LG3_572899434) were associated with SPAD in at least three out of the five environments, and on average each marker explained 7%-13% of the phenotypic variance (PV) measured as the difference in R-square of the model with the SNP and without the SNP. SNP markers  Note: soil plant analysis development (SPAD), spectral reflectance and canopy temperature were taken four to six times in a season during reproductive stage on hot days at solar noon. A SPAD reading > 50 indicates a dark-green color and high chlorophyll concentration, a reading < 40 indicates a yellow-green color and low chlorophyll concentration. Reproductive stem length, internode length and pod number were measured on three plants per plot at physiological maturity. The overall weather classification of environments 2015 and 2017 at Saskatoon was heat stress, and the remaining three environments condition was ambient (control) for pea production. A SPAD value is an index of light transmittance at 650 nm and 940 nm. Similarly, PRI is an index derived from narrow-band (531 and 571 nm) spectral reflectance.
Chlorophyll concentration, measured by a SPAD meter, was affected by genotype, environment and their interaction; and the variance component analysis showed that maximum variation (67.9%) among the GWAS panel was due to the genotype effect, and the broad sense heritability was 0.95. Overall, genotype chlorophyll concentration ranged from 26.6 to 57.6 SPAD values under heat stress, and 30.0 to 67.5 under control conditions (Table 3). On average, the heat stressed environments had 3% less SPAD value than the ambient environments. Six markers (Chr5LG3_150942510, Chr5LG3_446272814, Chr5LG3_449362407, Chr5LG3_566189589, Chr5LG3_569788697, and Chr5LG3_572899434) were associated with SPAD in at least three out of the five environments, and on average each marker explained 7%-13% of the phenotypic variance (PV) measured as the difference in R-square of the model with the SNP and without the SNP. SNP markers Chr5LG3_566189589 and Chr5LG3_449362407 were associated with SPAD in 4 and 5 environments explaining 13% and 11% of the PVs, respectively (Table 4). PRI was also significantly affected by genotype, environment and by the G x E interaction. Variance components showed most of the variation in PRI was due to environmental factors, and the broad sense heritability was the least (0.35) compared with the other traits ( Table 2). Two markers, Chr6LG2_469101917, and Chr7LG7_263964018 were significantly associated with PRI at three out of the five environments. Each of the two markers explained 9% of PV (Table 4).  Note: All markers presented here were significant in at least three of five environments for a given trait. In each SNP designation, Chr and LG indicate chromosome and linkage group and the number after the _ is the base pair position. For non-chromosomal SNPs, Sc refers to scaffold followed by the scaffold number. Each locus is represented by one SNP marker of the LD block. † R-square value is presented as the difference of R-square explained by the model with and without SNP.
For canopy temperature (CT), the GWAS accessions significantly varied due to both genotype (G) and environment (E) effects, but not by the G x E interaction ( Table 2). In general, under heat stress, the accessions' CT, measured four to six times in a season during reproductive stage on hot days at solar noon, ranged from 24.5 to 31.0 • C, whereas under ambient conditions, the CT ranged from 21.4 to 26.9 • C. This temperature difference indicated that CT is highly influenced by the environment effects with a relatively lower broad sense heritability of 0.57 (Table 2; Table 3; Figure 1). Two SNP markers (Chr4LG4_415994869 and Chr5LG3_309595819) were associated with CT in three of the five environments. The R-square value of the model with SNP ranged from 0.43 to 0.53, and each of the SNP markers explained 6% of PV.
Reproductive stem length was also affected by genotype and environment main effects and their interaction. The reproductive stem length under the stressed environments ranged from 13 to 99 cm, whereas under the control environments the range was from 14 to 117 cm, suggesting heat stress decreased the reproductive stem length. Analysis of variance components showed genotype and environment main effects respectively contributed to 63.4% and 7.6% of the variation in the GWAS panel. The broad sense heritability for reproductive stem length was 0.92. Seven SNP markers (Chr3LG5_18678117, Chr4LG4_29062302, Chr5LG3_566189271, Chr5LG3_572669963, Chr7LG7_20086906, Chr7LG7_23295474, and Chr7LG7_96157380) were associated with reproductive stem length in at least three of the five environments, and four of these SNPs were consistent in at least four of the five environments. SNP marker Chr4LG4_29062302 was found to be associated with the trait in all five environments with an average R-square of the model of 0.60. Overall, the R-square value of the model with SNP ranged up to 0.71 for reproductive stem length (Table 4).
Internode length was another trait significantly affected by genotype and environment main effects and their interaction. Under heat stress, the internode length ranged from 1.6 to 11.3 cm with a mean value of 11.0 cm, whereas under control conditions, the range was 1.9 to 14.9 cm with a mean value of 14.8 cm. Variance component analysis showed genotype and environment respectively contributed 43% and 4.8% of the variations to the GWAS panel. The broad sense heritability was 0.90. Six SNP markers (Chr4LG4_62461234, Chr4LG4_63111072, Chr4LG4_80759704, Chr5LG3_566189271, Chr6LG2_420562729, and Chr7LG7_197862543) were associated with internode length in at least three of the five environments. These markers were significantly associated with internode length in at least three of the five environments with the R-square value of the model with SNP ranged up to 0.63. SNP marker Chr5LG3_566189271 was identified in all five environments with an average R-square of 0.49.
Pod number was also significantly affected by genotype and environment main effects and their interaction. Variance component analysis showed genotype and environment, respectively, contributed 36.6% and 12.4% to the overall pod number variance in the GWAS panel. Compared with the three control environments, pod number under the heat stress environments decreased by 14.6%. The broad sense heritability in pod number was 0.88. Eight SNP markers (Chr2LG1_4359822, Chr2LG1_105547608, Chr2LG1_370541780, Chr2LG1_385949935, Chr2LG1_389336188, Chr2LG1_402022079, Chr3LG5_216337201, Chr5LG3_530537682, and Sc04062_32372) were associated with pod number in at least three of the five environments explaining 7% to 9% of PV, with an average R-square value of 21.9.
Manhattan plots showing the association of SNP markers with plant chlorophyll concentration and reproductive stem length in multiple trials, and the corresponding Q-Q plots are presented as examples from this research in Figures 2 and 3, respectively. The Q-Q plots represent the observed P values of each SNP marker against the expected P values. The Manhattan plots in Figure 2 showed the significant association of SNP markers on Chr 5 (LG3) with plant SPAD in each of the individual environments presented. The Manhattan plots in Figure 3 showed the significant association of SNP markers on multiple chromosomes with the reproductive stem length.

Overall Association of Phenotypic Traits
Principal component analysis (PCA) based on the correlation of traits revealed the overall traits association and the genotype response across the five environments ( Figure 4A,B). The first two PCs explained 61.9% of the total variability in the data. The loading plot illustrated traits association and how much each trait contributed to the PCs. The first PC was influenced mainly by SPAD, reproductive stem and internode lengths, whereas the second PC was influenced mainly by CT and Of all the MTAs that were observed in > 60% of the environments, the following markers had the greatest percent variation averaged over the selected environments for the respective traits: Chr5LG3_566189589 (13% PV) and Chr5LG3_449362407 (11% PV) for SPAD; Chr6LG2_469101917 and Chr7LG7_263964018 each with 9% PV for PRI; Chr4LG4_415994869 and Chr5LG3_309595819 each with 6% PV for CT; Chr3LG5_18678117 (6% PV), Chr5LG3_572669963 (5% PV), and Chr7LG7_96157380 (5% PV) for reproductive stem length; Chr4LG4_63111072 (6% PV), Chr5LG3_566189271 (7% PV) and Chr4LG4_62461234 (6% PV) for internode length; and seven markers, Chr2LG1_105547608, Chr2LG1_370541780, Chr2LG1_385949935, Chr2LG1_389336188, Chr3LG5_216337201, Chr5LG3_530537682, and Sc04062_32372 each with 9% PV for pod number (Table 4).
Forty-eight unique genes were identified within a 15 kb region of the selected 32 SNP markers and are considered as candidate genes. The candidate genes identified for various traits included those encoding for transcription factor, translation initiation factor, heat shock protein, ribosomal protein, protein kinase, transmembrane protein, and others as listed in Table 5. Two genes, Psat5g299080 and Psat5g299040, which encode the proteins kinesin-related protein 4-like and PPR containing plant-like protein (putative tetratricopeptide-like helical domain-containing protein), were identified as potential candidate genes associated with internode length, reproductive stem length and chlorophyll content (SPAD).

Overall Association of Phenotypic Traits
Principal component analysis (PCA) based on the correlation of traits revealed the overall traits association and the genotype response across the five environments ( Figure 4A,B). The first two PCs explained 61.9% of the total variability in the data. The loading plot illustrated traits association and how much each trait contributed to the PCs. The first PC was influenced mainly by SPAD, reproductive stem and internode lengths, whereas the second PC was influenced mainly by CT and pod number. SPAD positioned in an opposing direction (obtuse angle to straight line) to reproductive stem and internode lengths indicating a significant negative correlation between SPAD and the length measurements. Likewise, CT positioned in the opposite direction of pod number indicating their significant negative correlation. The hotter the canopy, the lower the pod number and thus seed yield ( Figure 4A). Score plots illustrated genotype placement (response) across the environments ( Figure 4B). The heat and or drought stressed environments (2015 Saskatoon and 2017 Saskatoon) positioned to the negative direction PC2 associating with high CT, whereas the control environments were associated greater pod number and SPAD value.

Discussion
As a cool season crop, pea is sensitive to heat stress which causes a significant yield loss. However, there exists substantial genetic variation among pea genotypes for heat tolerance [10,12,24,28]. A strategic assessment and use of available variation is essential for crop improvement through using allelic variation. With the availability of cost-effective, high-throughput SNP genotyping methods and genomic resources, GWAS has been an effective method for identifying genetic loci associated with traits of many crop species including legumes [29,30,36].
The present GWAS was undertaken to identify SNP markers associated with traits related with pea heat response using a panel of 135 genetically diverse pea accessions. The accessions were from breeding programs of major pea growing areas and, thus accounted genotypes with a wide range of heat sensitivity. Genotyping by sequencing (GBS) identified 16,877 good quality SNPs, of which 15,609 were distributed across seven chromosomes of pea and the remaining 1268 were nonchromosomal SNPs [30].
Linkage disequilibrium patterns of population structure and genetic relatedness information are important for association mapping to minimize the number of false positive associations [41], thus the LD of the 135 GWAS members was analyzed by chromosome, and the LD decay estimates of the 7 chromosomes ranged from 0.03 to 0.18 Mb [30]. Based on genetic relatedness the 16,877 SNPs in the GWAS panel were clustered into nine groups [30]. Similarly, Diapari et al. [37] clustered another 94 pea accessions into eight groups, and Siol et al. [42] grouped 917 Pisum accessions into 16 groups. The above groupings indicated the extent of genetic variability among pea accessions. The clustering did not necessarily correspond solely with the geographic origins of the individuals, but depended on additional factors of variability such as the objectives in different breeding programs [30].
In the present GWAS, we evaluated ten heat stress-responsive traits. The first six were: chlorophyll concentration by SPAD, PRI, CT, reproductive stem length, internode length, and pod

Discussion
As a cool season crop, pea is sensitive to heat stress which causes a significant yield loss. However, there exists substantial genetic variation among pea genotypes for heat tolerance [10,12,24,28]. A strategic assessment and use of available variation is essential for crop improvement through using allelic variation. With the availability of cost-effective, high-throughput SNP genotyping methods and genomic resources, GWAS has been an effective method for identifying genetic loci associated with traits of many crop species including legumes [29,30,36].
The present GWAS was undertaken to identify SNP markers associated with traits related with pea heat response using a panel of 135 genetically diverse pea accessions. The accessions were from breeding programs of major pea growing areas and, thus accounted genotypes with a wide range of heat sensitivity. Genotyping by sequencing (GBS) identified 16,877 good quality SNPs, of which 15,609 were distributed across seven chromosomes of pea and the remaining 1268 were non-chromosomal SNPs [30].
Linkage disequilibrium patterns of population structure and genetic relatedness information are important for association mapping to minimize the number of false positive associations [41], thus the LD of the 135 GWAS members was analyzed by chromosome, and the LD decay estimates of the 7 chromosomes ranged from 0.03 to 0.18 Mb [30]. Based on genetic relatedness the 16,877 SNPs in the GWAS panel were clustered into nine groups [30]. Similarly, Diapari et al. [37] clustered another 94 pea accessions into eight groups, and Siol et al. [42] grouped 917 Pisum accessions into 16 groups. The above groupings indicated the extent of genetic variability among pea accessions. The clustering did not necessarily correspond solely with the geographic origins of the individuals, but depended on additional factors of variability such as the objectives in different breeding programs [30].
In the present GWAS, we evaluated ten heat stress-responsive traits. The first six were: chlorophyll concentration by SPAD, PRI, CT, reproductive stem length, internode length, and pod number. The other four were: plant height, lodging, pod to node ratio, and water band index (WBI). From the latter four traits, five SNP markers on Chr 1 (LG6), Chr 2 (LG1), Chr 3 (LG5), Chr 5 (LG3) for lodging, and four SNP markers on Chr5 (LG3) for plant height were previously reported by Gali et al. [30], and no marker was detected to be significant in at least three of the five environments for pod to node ratio and WBI. As such, in the current paper we focused on the first six traits for phenotypic variation in the 135 pea accessions across five environments.
The five environments were grouped into ambient (three environments) and heat and or drought stress (two environments) conditions based on weather data and threshold temperature for heat stress in the field [5]. All traits had a wide range of phenotypic variation within each environment and stress level, which is essential for dissecting complex traits through association mapping. Overall, we identified 32 MTAs for six traits that have physiological and agronomic importance and are involved in pea heat response. A marker identified for a significant association with a given trait would be more reliable if the same marker is found in multiple environments [30]. Therefore, for the six traits we investigated, the SNP markers deemed significant were consistent in at least three environments, and these markers could potentially be used for marker-assisted selection of these traits in the effort of improving pea for heat tolerance.
In this study, the SPAD value was used to estimate chlorophyll concentration, a major component of chloroplasts, and can be used as a factor to determine crop adaptation to environmental stresses by retention of greenness [10,13,43]. Regression analysis on wheat reported that under heat stress, the SPAD value was associated with plasma and thylakoid membrane damage [44], which hinders light absorbing efficiency of photosystems (PSI and PSII), and hence reduced photosynthetic capacity ultimately leading to crop yield loss [11]. Understanding of the genetic bases that govern chlorophyll concentration may contribute to enhancing photosynthetic efficiency and thus minimize yield loss due to stressful environments.
We identified six MTAs that were related to SPAD value in repeated tests. All of the MTAs identified for SPAD were from Chr 5 (LG3). Bell et al. [45] reported that pea chlorophyll degradation under stress conditions is governed by the SGRL protein, a distinct class of the SGR gene which is induced by environmental stresses. The genomic location of SGRL was reported to be on LG3 which supports our result where all of the SPAD markers also reside on Chr 5 (LG3). The SGRL gene sequence (https://www.ncbi.nlm.nih.gov/nuccore/LN810021) location in the pea genome assembly spanned between the base pair positions Chr5LG3_151800929 and Chr5LG3_151804253, and is within close proximity (858 Kbp) of the SPAD associated marker ChrLG3_150942510. Using GWAS on soybean, Dhanapal et al. [29] identified 52 SNP markers associated with chlorophyll content.
Similarly, two loci were identified for association with PRI. One of these loci is on Chr 6 (LG2) and the second is on Chr 7 (LG7), and in both cases the markers were consistent in three environments. There are only a few reports that have used GWAS to identify markers associated with vegetation indices, namely, in soybean and wheat. In soybean, Herritt et al. [25] identified 31 SNPs linked with PRI, and on wheat, Gizaw et al. [34] reported the presence of markers associated with normalized chlorophyll-pigment ratio index (NCPI), and normalized difference vegetation index (NDVI). However, use of GWAS and vegetation indices has been lacking in cool season pulse crops. To the best of our knowledge, our report is the first to apply VIs in pea GWAS. The PRI is increasingly used as a predictor of crop photosynthetic efficiency which responds to environmental variables [19]. PRI is associated with photosynthetic protective mechanisms by dissipation of excess energy such as in the operation of xanthophyll cycle during stress. Violaxanthin de-epoxidase VDE is among the genes known to be involved in excess energy dissipation in the xanthophyll cycle [46]. Two MTAs, one each on Chr 5 (LG3) and Chr 4 (LG4), were detected for CT, a trait consistently used as an indicator of stress mainly of drought and heat stresses [12,14]. Generally, cooler canopy is associated with heat avoidance, and is an indicator of a healthy canopy with an optimal physiological state [12]. Again, to the best of our knowledge, no previous study exists on pea or other cool season legume crops that has reported genomic regions associated with CT. In a study using 24 pea cultivars across six environments, Tafesse [14] reported that leaf surface wax concentration is positively correlated with water band index, a proxy for leaf water retention, and contributes to a cooler canopy. WAX2 is among the genes controlling wax biosynthesis in Arabidopsis [47], and glossy13 is another gene with similar role reported in maize [48]. Lodging contributes to canopy heating in pea, and upright and semileafless cultivars with the afila gene have cooler CT [12,49]. Tar'an et al. [50] identified major loci associated with lodging resistance in pea on LG III, and one of the markers we identified for CT is also on LG III, suggesting genes controlling lodging also control CT.
We identified seven MTAs associated with reproductive stem length on chromosomes 3 (LG5), 4 (LG4), 5 (LG3), and 7 (LG7); and six MTAs associated with internode length on chromosomes 4 (LG4), 5 (LG3), 6 (LG2) and 7 (LG7). The markers associated with these two traits mostly were positioned on the same linkage groups, and a SNP marker Chr5LG3_566189271 was associated with both the traits. Using the current GWAS panel, Gali et al. [30] identified four MTAs associated with plant height that were on same linkage group as that of reproductive stem length and internode length. The SNP marker Chr5LG3_566189271 reported for plant height [34] was also associated with internode length in the current study. Both reproductive stem and internode lengths were significantly reduced by heat stress [12]. A cultivar's genetics affects internode length, and in pea the Le gene controls internode length [25], which directly affects reproductive stem length and plant height via its influence on gibberellic acid function on growth and determinacy/indeterminacy [25,51,52]. Using two pea recombinant inbred populations, Weeden [22] identified a major QTL on LG3 for a longer internode (Le), and a second QTL on LG4 for the recessive allele which caused plants to have shorter internodes.
We identified nine loci associated with pod number, of which six were on Chr 2 (LG1), one each on Chr 3 (LG5) and 5 (LG3), and one on a non-chromosomal scaffold. Plant pod number is the number of flower-bearing nodes multiplied by the average number of flowers per node. Previously, Jiang et al. [28] identified two unmapped QTLs for pod number using 92 diverse accessions. Also, Huang et al. [24] identified two QTLs for pod number based on a bi-parental mapping population on Chr 5 (LG3). The greater number of loci identified in this study was likely due to the use of a GWAS panel which represented a broad range of diversity in pod number ranging from 3 to 19 pods per plant, and where most of this variation is contributed by genetic factors. Benlloch [53] indicated that flower number per plant, which directly determines pod number, is controlled by two genes, Fn and Fna, and a single mutation of these genes increases flower number per plant. Pod number is a major yield component that has a strong correlation with seed yield, and is most affected by heat stress [12,23,24]. The reduction in pod number and yield was likely from heat stress-induced abortion of flower buds, flowers, and pods [4,23]. Pod set relies on pollen and stigma functioning optimally, both of which are very sensitive to heat stress [54].
In conclusion, in this GWAS we identified 32 MTAs and 48 candidate genes for traits associated with pea heat response. These results are expected to enhance the understanding of genetic loci controlling these traits. The identified candidate genes are involved in various biological functions and require further functional validation. The detected MTAs and candidate genes should be useful for marker-assisted selection for heat tolerant pea varieties.

Plant Materials
A panel of 135 diverse field pea accessions, as described by Gali et al. [30], were grown for two years (2016-2017) at two Rosthern (52 • 19 were from Australian pulse breeding programs, 77 were from eastern and western European countries, the Russian Federation and the UK, 15 were from the USA, 17 were from Canada (mostly from the Crop Development Centre, University of Saskatchewan), five were from Ethiopia, and two were from India. Thus, the accessions represented the major pea growing areas of the world. The accessions were commercial cultivars released over the past 50 years for local production, and were able to flower and mature under the five environments tested [30].

The Field Trials and Weather Conditions
The experimental design at each environment was a randomized complete block with two replications. Plot size was 1.37 m width × 3.66 m length, and the recommended seeding rate (100 seeds m −2 , targeting 80-85 plants m −2 on 0.25 m row spacing) was used. Weed control was achieved by management practices used in pea production in Saskatchewan as described by Tafesse et al. [12].
Weather data for 2015 Saskatoon starting from June 11 to the end of the growing season, 2016 Rosthern starting from June 21 to the end of the growing season, and 2016 Saskatoon starting from July 21 to the end of the season were collected from weather stations (Coastal Environmental Systems, Seattle, WA, USA) established at each site. Weather data of 2017 and the remaining 2015 and 2016 were obtained from Environment Canada database (https://climate.weather.gc.ca) recorded by the nearest stations to the trial sites. For Saskatoon, data from central Saskatoon station, and for Rosthern the mean of data from Saskatoon international airport and Prince Albert stations were used. The daily maximum air temperatures, amount of precipitation and number of days when air temperature exceeded 28 • C during the growing seasons were used to determine the degree of stress in each environment at different growth stages. The categorization of growth stages into vegetative (germination to end of vegetative growth) and reproductive (beginning of flowering to maturity) was conducted using the phenology data reported by Gali et al. [30]. Based on the weather data, 2015 and 2017 Saskatoon had heat and drought stress conditions and the remaining three environments were generally ambient and considered control environments (Table 1).

Phenotypic Measurements
Chlorophyll concentration was estimated non-destructively using a SPAD502Plus chlorophyll meter (Konica Minolta Sensing Americas Inc., USA). The SPAD value is a unitless index, calculated as the ratio of the intensity of light transmittance at red (650 nm) to infrared (940 nm) and gives a value that corresponds to the relative amount of chlorophyll present in the leaf. Hereafter, the chlorophyll concentration estimated by SPAD meter is referred to as 'SPAD'. The SPAD readings were taken four to six times each season, and for each measurement day the mean SPAD value was calculated by the instrument from three readings taken from three plants per plot on fully expanded stipules at the second or third node counting down from the apex of a main stem.
Similarly, spectral measurement was conducted repeatedly on leaf stipules using a portable spectroradiometer PSR-1100F (Spectral Evolution Inc, Lawrence, MA, USA). This device enabled hyperspectral readings with a range of 320-1,126 nm, and 1.6 nm sampling interval, and a total of 512 discrete narrow bands. PRI was calculated from the reflectance data according to Gamon et al. [19] as: PRI = (R 531 − R 570 )/(R 531 + R 570 ) (1) where R is reflectance percentage and 531 and 570 are the wavelength bands in nm along the light spectrum. The PRI is used as a proxy for the xanthophyll cycle, a photosynthetic protective cycle that operates more during stress [19].
Canopy temperature (CT) was measured four to six times in each location in a season using a hand held infrared thermometer (Model 6110.4ZL, Everest Interscience Inc, Tucson, AZ, USA) as described by Tafesse et al. [12]. Measurements of SPAD, spectral reflectance, and CT were carried out repeatedly (four to six times in a season) during the reproductive stage, at solar noon on relatively hot days when air temperature is greater than 25 • C, and the mean value was used for analysis.
The other measurements taken at physiological maturity were: reproductive stem length (vine length from first flowering node to the tip of the main stem); internode length (determined as the ratio of reproductive stem length to reproductive node numbers); and pod number per plant (total pods counting all pods with at least one seed on the main stem). For these, each measured variable was the mean of three plants per plot sampled at random and lengths were measured in cm.

Phenotype Data Analysis
Before employing analysis of variance (ANOVA), homogeneity of variances and normality of residuals were tested using checked using Levene and Shapiro-Wilk tests, respectively [55,56]. Variance components of genotype, environment, the G × E interaction, block within environment, and the residual were analyzed using the generalized linear model (GLM) and by considering all factors as random effects. Broad sense heritability (H 2 ) was calculated as: where σg 2 is the genetic variance, and σge 2 is the variance of genotype and environment [57]. Over environments, combined ANOVA on SPAD, PRI, CT, reproductive stem length, internode length, and pod number was carried out using the Mixed procedure of SAS (Version 9.4, SAS Institute). Genotype, environment and G x E interaction were considered as fixed, and blocks as random factors. Principal component analysis (PCA) was performed with the multivariate function of Minitab (Version 19, Minitab LLC, USA) using means of traits to infer overall association among traits and genotype for the five environments. Based on significant eigenvalue (> 1), the first two principal components (PC) were selected for the minimum number of PCs to explain the greatest total variation in the data set.

Association Mapping
Genotyping of the 135 GWAS panel was performed by genotyping-by-sequencing (GBS, [58]), and 16,877 SNPs were reported based on a minimum read depth of five and minimum allele frequency of 0.05 [30]. The reported SNPs were used for association analysis using GAPIT (Genome Association and Prediction Integrated Tool-R package [30]) software. Association analysis for each trait was conducted using the mixed linear model (MLM). For MLM analysis, Q values were generated from structure analysis [59] and K (kinship coefficient matrix) values calculated by GAPIT and identity-by-state (IBS) methods were used. Principal co-ordinate values were used as co-variates. Although the result is not presented here, the model output of MLM was compared with the Super MLM model and the markers identified in both methods were mostly similar. The quantile-quantile (Q-Q) plots of each trait were drawn using the observed and expected log 10 P values. Marker-trait associations were selected based on P value (P ≤ 0.001) and repeated occurrence of the association in at least three of the five trials. The genes within 15 kb of the identified markers are reported as the candidate genes. The pea genome sequence reported by Kreplak et al. [40] was used for identification of candidate genes.