Marker-trait association analysis for drought tolerance in smooth bromegrass

Abstract

Background

Little information is available on the application of marker-trait association (MTA) analysis for traits related to drought tolerance in smooth bromegrass. The objectives of this study were to identify marker loci associated with important agronomic traits and drought tolerance indices as well as fining stable associations in a diverse panel of polycross derived genotypes of smooth bromegrass. Phenotypic evaluations were performed at two irrigation regimes (normal and deficit irrigation) during 2 years; and association analysis was done with 626 SRAP markers.

Results

The results of population structure analysis identified three main subpopulations possessing significant genetic differences. Under normal irrigation, 68 and 57 marker-trait associations were identified using general linear model (GLM) and mixed linear mode1 (MLM), respectively. While under deficit irrigation, 61 and 54 markers were associated with the genes controlling the studied traits, based on these two models, respectively. Some of the markers were associated with more than one trait. It was revealed that markers Me1/Em5–11, Me1/Em3–15, and Me5/Em4–7 were consistently linked with drought-tolerance indices.

Conclusion

Following marker validation, the MTAs reported in this panel could be useful tools to initiate marker-assisted selection (MAS) and targeted trait introgression of smooth bromegrass under normal and deficit irrigation regimes, and possibly fine mapping and cloning of the underlying genes and QTLs.

Background

Smooth bromegrass (Bromus inermis Leyss.) is a long-lived and cool-season grass species [1, 2] which is adapted to dry regions and grown mostly for hay, pasture and soil conservation [1].

Drought is the major abiotic constraint affecting growth and productivity of crops worldwide [3, 4]. Development of suitable cultivars with more drought tolerance is crucial for enhancing sustainable production of plants and provides a strategy for alleviation of the effects of climate change [5]. However, drought tolerance is a complex and quantitative trait, including several metabolic pathways. A promising strategy to facilitate selection for drought tolerance is to identify and select for genetic markers linked to the traits related to drought tolerance (marker assisted selection; MAS). The main prerequisite for MAS is the availability of markers that are tightly linked to genes or quantitative trait loci (QTLs) which can be used to select for traits that are difficult to measure or dependent on the developmental stage [6, 7].

The application of molecular markers allows locating the genes of interest in the genome, thus avoiding the time and the space needed in breeding programs [8]. Sequence-related amplified polymorphism (SRAP) is a PCR-based molecular marker technique which can be used for a variety of purposes, including genetic diversity analysis, map construction, gene tagging, genomic and cDNA fingerprinting, and map based cloning [9, 10]. Moreover, it is an advanced molecular marker for genetic research in grass and forage species, which uses from two primers with an arbitrary nucleotide sequence and therefore can detect nucleotide sequence polymorphisms [11]. Analysis of SRAP data has frequently been used for the construction of linkage maps [12, 13] and identification of QTLs [14, 15]. However, the utilization of SRAP for grass and forage researches such as association analysis in smooth bromegrass is rare.

The first step towards MAS, as an important tool for accelerating varietal improvement and rate of genetic gain, is to dissect marker-trait associations (MTAs) [16]. Two approaches which can be used for dissection of quantitative complex traits are association mapping and linkage analysis [17, 18]. However, QTL mapping has a low resolution and requires a lot of time and resources [19, 20]. Genome-wide association studies (GWAS) or association mapping have recently become a popular alternative to QTL mapping for identifying MTAs in plant populations [21, 22]. Compared to linkage mapping, GWAS overcomes several of the drawbacks of QTL mapping. It offers higher mapping resolution, is less time consuming and requires fewer resources, and evaluates a wide range of alleles rapidly [23, 24].

In order to avoid identifying spurious associations between markers and traits and also to avoid of both types of error (types I and II) in association analysis, it is necessary to evaluate the population structure and use from a mixed-model approach [6, 25, 26]. Two models of general linear model (GLM) and mixed linear model (MLM) are used for association analysis. In MLM, both the kinship matrix (K) and population structure (Q) are merged, whereas in the GLM, only population structure information is utilized as a covariate [6].

In recent years, the application of association analysis in forage grasses is discussed in several reports. In perennial ryegrass (Lolium perenne) the application of association mapping for some traits such as flowering time, leaf length, carbohydrate content, submergence tolerance, salinity and drought tolerance has been evaluated [27,28,29,30,31]. In tall fescue (Festuca arundinacea), SSR loci related to agronomic traits [32] and heat-tolerance-related traits [33] have been detected. Kempf et al. (2017) applied SRAP and SSR markers in marker–trait association analysis for agronomic and compositional traits of sainfoin (Onobrychis viciifolia) [7]. In orchardgrass, studies have been carried out for detection of the loci related to drought tolerance, seed yield, forage yield [34], rust resistance [35], and heading date [36]. Such studies have demonstrated that GWAS is an efficient method for identifying genomic regions associated with complex quantitative traits. However, in smooth bromegrass the application of association mapping in identifying links between genes or markers with complex quantitative traits such as drought tolerance is still in its infancy. This study was conducted to: i) identify genetic loci associated with the key agronomic traits and drought tolerance, under normal and water deficit conditions using SRAP markers; and ii) discover stable marker loci linked to the agronomic and drought tolerance related traits.

Results

Phenotyping

Significant differences were observed between irrigation regimes for most of the measured traits. Except for flag leaf length (FLL), flag leaf width (FLW), panicle length (PL), and winter growth vigor (WGV), the magnitude of mean performance was significantly decreased for all of the evaluated traits, under water-deficit condition. Dry matter yield (DMY) was decreased by deficit irrigation 42% on average (Table 1).

Table 1 Mean performance, phenotypic coefficient of variation (PCV), genotypic coefficient of variation (GCV), and broad- sense heritability (h²_b) of traits recorded under normal and water deficit conditions in smooth bromegrass genotypes

Phenotypic coefficient of variation (PCV) showed a range from 4.89% for plant height (PH) to 23.77% for number of stems per plant (NS) under normal irrigation and from 6.89% for days to anthesis (DA) to 27.42% for NS under deficit irrigation (Table 1). Genetic coefficient of variation (GCV) ranged from 6.08% for PH to 27.32% for NS under normal irrigation and from 7.07% for DA to 30.87% for NS under deficit irrigation. Except for crown diameter (CD), the values of genetic variation under deficit irrigation were higher than the ones for normal irrigation (Table 1).

The estimates of broad-sense heritability for each irrigation regime are given in Table 1. Moderate to high values of heritability estimates were observed for all of the evaluated traits, under both irrigation regimes. The range of this parameter was from 61.33% for WGV to 94.05% for DA under normal irrigation and from 66.01% for DMY to 95.13% for days to panicle emergence (DPE) under deficit irrigation regime. For all traits, estimates of heritability were higher under normal irrigation than deficit irrigation (Table 1).

Genotyping

In total, 959 bands were created from 30 SRAP primer combinations, of which 626 bands were polymorphic (Table 2). The range of total number of bands scored per primer combination was from 13 (Me3/Em3) to 28 (Me4/Em2), with an average of 21 bands. The percentage of polymorphic bands ranged from 55.17% (Me4/Em3) to 85.19% (Me3/Em4) with an average of 69.85%. The relative informativeness of each marker can be evaluated based on its PIC value. In the present study, PIC value ranged from 0.35 (Me1/Em5) to 0.50 (Me2/Em4, Me3/Em4, and Me5/Em3), with the average of 0.45. The highest and lowest MI values were 12.60 (Me4/Em2) and 5.85 (Me3/Em3), respectively. Markers Me4/Em2 and Me3/Em1 showed the highest and lowest RP values, respectively (Table 2).

Table 2 Information and diversity statistics for sequence related amplified polymorphism (SRAP) markers used for association analysis in smooth bromegrass

Population structure and association analysis

The maximum likelihood and DK were used to calculate the optimum number of subpopulations (K). The maximum value of DK obtained at K = 3, suggested that there were three subpopulations in the smooth bromegrass panel (Table 3; Figs. 1 and 2). Each of these three subpopulations had its own characteristics. Subpopulation 1 contained seven genotypes of G30, G31, G32, G33, G34, G35, and G36 (Fig. 2). All of these genotypes were belonged to Hungary, and had higher persistence than other genotypes. Subpopulation 2 included genotypes G11, G15, G21, G26, and G28, which all of them were belonged to Hungary except for G15 (Iranian genotype). Genotypes of this subpopulation were early flowering and had lower productivity than other genotypes. The remaining genotypes were located in the third subpopulation (Fig. 2). This subpopulation showed late flowering and had higher productivity than other ones. All of the genotypes of this subpopulation were Iranian except for G12, G13, and G25 (Hungarian genotypes). As observed, structure analysis was able to separate genotypes based on geographical or ecological data.

Table 3 Calculated statistics to detect optimum number of subpopulations (K) in structure analysis of smooth bromegrass genotypes using the STRUCTURE program

Fig. 1

Population structure analysis in a diverse germplasm of smooth bromegrass (The optimum number of subpopulations was determined using Δk in the Bayesian clustering method)

Fig. 2

The results of genetic association analysis of smooth bromegrass genotypes performed by STRUCTURE software version 2.3.4. The membership coefficients of genotypes are given on the y-axis. Each of the three identified subpopulations is shown by a different color

Association analysis between SRAP markers and the phenotypic mean of traits (over years and irrigation regimes) was separately conducted based on both GLM and MLM models. Under normal irrigation, based on the GLM model (P values < 0.01 and a cut-off value of 0.05 for the FDR) 68 SRAP markers showed significant associations with means of the studied traits, at the 0.01 probability level (Table S1). The percentage of phenotypic variation (coefficient of determination, R²) of an individual trait explained ranged from 11.75 to 31.32% (Table S1). Under deficit irrigation, 61 markers had significant associations with the studied traits, at the 0.01 probability level (Table S1). The percentage of phenotypic variation (R²) of a trait explained varied from 9.89 to 26.28% (Table S1). However, in the MLM model, kinship or relatedness matrix was considered as a factor, and the number of significant markers decreased as compared to GLM model. So that, under normal irrigation 57 markers and under deficit irrigation 54 markers indicated significant associations at the 0.01 probability level (Table 4). In this model, the percentage of phenotypic variation, under normal irrigation ranged from 7.71 to 20.89% and under deficit irrigation varied from 6.76 to 17.61% (Table 4). Moreover, association analysis was also done for drought tolerance and susceptibility indices. Results revealed that 19 and 20 markers showed significant associations with the calculated indices based on GLM and MLM model, respectively (Tables 5 and S2).

Table 4 Association of SRAP markers with phenological, morphological, and agronomic traits of smooth bromegrass genotypes under normal and water deficit conditions based on mixed linear model (MLM)

Table 5 Association of SRAP markers with drought tolerance and susceptibility indices of smooth bromegrass genotypes based on mixed linear model (MLM)

Based on the results of GLM and MLM models, some markers were associated with more than one trait at the same time. For example, under normal and deficit irrigation regimes, marker Me1/Em6–7 showed simultaneously significant associations with DA, PH, FLL, and FLW, based on both GLM and MLM models. Marker Me2/Em5–21 had concurrently significant associations with DA, PH, and FLL, under both irrigation regimes and in both GLM and MLM models (Tables 4 and S1). In addition, under normal irrigation, marker Me2/Em1–14 showed significant associations with DPE, DA, PH, FLL, and FLW, based on GLM model; and showed significant associations with the traits of DPE, DA, PH, and FLL based on MLM model. However, under deficit irrigation this marker had significant associations with DA and PH based on both models (Tables 4 and S1). Based on MLM model, marker Me1/Em5–11 showed associations with DMY and WGV under both normal and deficit irrigation regimes (Table 4). Moreover, markers Me1/Em5–11 and Me1/Em3–15 showed significant associations with MP, GMP, and STI, based on both GLM and MLM models (Tables 5 and S2).

To assess stable associations, association analysis was conducted in each irrigation regime, separately. In total, 30 and 21 MTAs showed adequately stable expression across irrigation regimes, based on GLM and MLM models, respectively. For instance, in GLM model, markers Me1/Em6–7, Me2/Em1–14, Me2/Em2–13, Me4/Em6–22, Me2/Em5–21, Me2/Em3–19, Me2/Em1–12, and Me1/Em6–16 showed significant and stable associations in both irrigation regimes with DA. Similarly, markers Me1/Em5–11, Me1/Em5–24, Me2/Em6–16, Me5/Em2–20, and Me1/Em2–19 displayed significant and constant associations in both irrigation regimes with WGV (Table S1). In MLM model, markers Me1/Em6–7, Me2/Em1–14, Me2/Em2–13, Me4/Em6–22, and Me2/Em5–21 had significant and stable associations in both irrigation regimes with DA. Also, markers Me4/Em4–14, Me4/Em6–2, Me5/Em2–1, Me1/Em2–21, and Me1/Em1–13 were constantly associated with FLW in both moisture conditions. In the same way, marker Me1/Em5–11 was associated with DMY (Table 4).

Discussion

Significant genetic variations among genotypes in terms of all of the evaluated traits demonstrate the difference in genes controlling yield, its components, and drought-tolerance. Moreover, the non-static performance of genotypes in two irrigation regimes emphasizes the importance of marker-trait association analysis in the two moisture environments, separately.

Most of the evaluated traits were significantly affected by water deficit due to decreased water potential of the soil and decline in net assimilation and photosynthesis of leaves [37, 38]. Similar results were reported in other researches [39, 40]. Wide genetic variation observed for all of the evaluated traits is promising genetic progress for these genotypes. Moreover, higher estimates for PCV and GCV under the deficit irrigation regime compared with normal irrigation indicate that water deficit have increased genetic variation for most of the studied traits and therefore, selection under deficit irrigation would be more effective. Our findings in this case are consistent with the reports of Abtahi et al. [41] and Saeidnia et al. [42]. However, some researchers stated that the genetic advance through selection is higher under normal irrigation [39, 43]. Moderate to high heritability estimates observed for all of the evaluated traits indicates that the improvement of these traits would be possible through selection and also emphasizes that detecting of marker–trait associations is possible for these traits [44].

The high percentage of polymorphism indicated that SRAP markers used in the present research could be used as powerful tools for discriminating of smooth bromegrass genotypes. Results of this study also revealed that the primer combination Me4/Em2 with the highest polymorphism percentage and high values of PIC, MI, and RP indices is informative and powerful enough for identification and discrimination of smooth bromegrass genotypes.

Population structure analysis identified three groups of genotypes in the studied panel of smooth bromegrass. As expected, structure analysis was able to separate genotypes based on their origin. Based on the results of association analysis of different traits under normal and deficit irrigations, the number of significant MTAs was lower in the MLM than GLM model. SRAP markers identified based on the results of MLM model can be considered as the most interesting candidates for future studies using MAS. It is stated that the combination of Q and K matrices strongly reduces the coefficient of determination and likely provides the best correction for population structure [45].

Marker-trait associations were mostly different in normal and deficit irrigation regimes. The results showed that a greater number of genes were probably involved in controlling traits at deficit irrigation regime than normal one. The percentage of variation which is explained by identified associations was low (7.71–20.89% under normal irrigation and 6.76–17.61% under water-deficit irrigation). This low R² value for each trait may be attributed to the role of many minor genes controlling the trait, outcrossing nature of smooth bromegrass, markers exhibiting minor quantitative effect, rare alleles, and complex allelic interactions [46, 47]. These results are in agreement with the findings of Lou et al. [32] and Sun et al. [33] in tall fescue.

Based on the results of GLM and MLM models, some markers had simultaneously significant associations with more than one trait. These markers may be effectively used to improve several traits, concurrently [33, 34]. Multi-association between different traits could be attributed to the co-expression mediated by expression of quantitative trait loci or e-QTLs [48]. For instance, marker Me1/Em6–7 showed simultaneously significant associations with DA, PH, FLL, and FLW, under both irrigation regimes and based on both GLM and MLM models. Similarly, marker Me2/Em5–21 concurrently showed significant associations with DA, PH, and FLL. In addition, under normal irrigation, marker Me2/Em1–14 showed significant associations with DPE, DA, PH, FLL, and FLW, based on GLM model; and also showed significant associations with DPE, DA, PH, and FLL based on MLM model. At deficit irrigation, this marker had significant associations with DA and PH based on both models. These simultaneous associations of markers with multiple traits may be attributed to pleiotropic effects or to several tightly linked genes that affect multiple traits [33, 49].

Determination of the genetic basis of drought tolerance requires correlating the occurrence of molecular markers with phenotypic scores for prediction of DNA genomic regions which involves effective factors on the response of plants [50]. Marker–trait association analysis identified 19 and 20 loci related with drought tolerance and susceptibility indices based on GLM and MLM models, respectively. Among these, SRAP markers Me1/Em5–11 and Me1/Em3–15 showed significant associations with MP, GMP, and STI, based on both GLM and MLM models. Moreover, in both models, markers Me5/Em5–9 and Me5/Em3–10 showed a significant association with DSI. If the effectiveness of these regions in the genetic control of drought tolerance is confirmed, these markers could be potentially used for the improvement of drought tolerance in smooth bromegrass.

Most of the MTAs were different under normal and deficit irrigations, indicating that the environmental factors have affected these associations [51]. These results showed that different genes may be effective on the same trait in different environments [52] or there might be a change in the expression level of the same gene between the two environments [48]. In the present study, 21 markers showed stable associations with different traits under both irrigation regimes. Diapari et al. [53] stated that associated markers which were detected in two or more different environments are more reliable than those present in only one environment.

Conclusion

In conclusion, the efficiency of association analysis approach as a powerful tool for identifying and detecting genes and markers linked to complex traits of agricultural and economic importance was illustrated. Satisfactory levels of polymorphism were observed for the studied traits in the polycrossed population. Three subpopulations were identified in smooth bromegrass genotypes; and 90 significant MTAs were detected using GLM and MLM models, under contrasting water conditions. Among these, three MTAs were identified for drought tolerance. Moreover, it was demonstrated that SRAP markers can be used in the future breeding programs to enhance drought tolerance of smooth bromegrass. Some SRAP markers were associated with the key agronomic traits of this species. Environmental specificity of MTAs shows that genotype × environment interactions are effective on association analysis; nevertheless, 30 and 21 MTAs showed significantly stable expression across two irrigation regimes based on GLM and MLM models, respectively. The markers identified in the present study are useful genomic resources for MAS in the future breeding programs of smooth bromegrass.

Methods

Plant materials and field experiment

A replicated nursery of smooth bromegrass (containing 1000 samples) was established at Isfahan University of Technology Research Farm in 2006. Some of these plant materials were provided by the Hungarian Institute for Agrobotany (HIFA), Tapioszele, Hungary. The Iranian ecotypes were natural ecotypes collected from wide geographical areas of Iran by a team which are specialist in the field of grass species. These samples were evaluated at first in the research center of Fozveh, Isfahan, Iran. After formal identification and verification of them, these ecotypes were used for research projects of Isfahan University of Technology (IUT).

Genotype panel used for the present study consisted of 216 clones randomly selected from a large nursery, comprised of 1800 single spaced-plant polycrossed progenies resulting from 36 parental ecotypes of smooth bromegrass (Bromus inermis). These 36 genotypes were randomly selected from polycross progenies of a set of 25 parental genotypes (Table 6). The parental genotypes were randomly chosen from the above mentioned nursery. For the present study, 216 clones were propagated in a greenhouse during the winter of 2012 and then were space-planted (50-cm grid) in the field according to a RCBD with six replications. Genotypes were evaluated under two levels of irrigation including a normal and a water deficit condition for 2014–2015. Under normal and water deficit conditions, irrigation was done when 50 and 90% of the total available soil water was exhausted from the root zone, respectively [54].

Table 6 Information on the genetic materials used in this study

This experiment was carried out in the field at Research Farm of Isfahan University of Technology, situated in Lavark, Najaf-Abad, Isfahan, Iran (32° 30′ N, 51° 20′ E, 1630 m amsl) during 2 years. This region has a mean annual temperature of 14.6 °C and mean annual precipitation of 141 mm, generally without rain during the summer, making supplemental irrigation necessary for growing crops during this period.

Phenotyping

During the year that plants were established (2013) no data was recorded. Days to panicle emergence (DPE), days to anthesis (DA), plant height (PH), and winter growth vigor (WGV) were measured as mentioned in our previous studies [39, 55]. Flag leaf length (FLL), flag leaf width (FLW), panicle length (PL), and number of stems per plant (NS) were measured at the pollination stage. After the complete flowering (about early summer), the produced forage of each genotype was harvested by cutting the grass from 5 cm above the ground in late spring, late summer, and late autumn. The harvested materials were then dried at 75 °C for 48 h and weighed to obtain dry matter yield (DMY). The average forage weight (g per plant) from the three cuts was used for analysis. The width of plant basal cover remained after the first harvest was considered as crown diameter (CD).

Five selection indices including tolerance index (TOL) [56], mean productivity (MP) [56], geometric mean productivity (GMP) [57], drought susceptibility index (DSI) [58], and stress tolerance index (STI) [57] were calculated based on the dry matter yield under normal and water deficit irrigations, according to the following formulations:

$$ \mathrm{TOL}={\mathrm{Y}}_{\mathrm{pi}}-{\mathrm{Y}}_{\mathrm{si}} $$

$$ \mathrm{MP}=\left({\mathrm{Y}}_{\mathrm{pi}}+{\mathrm{Y}}_{\mathrm{si}}\right)/2 $$

$$ \mathrm{GMP}={\left({\mathrm{Y}}_{\mathrm{pi}}\times {\mathrm{Y}}_{\mathrm{si}}\right)}^{0.5} $$

$$ \mathrm{DSI}=\left[1-\left({\mathrm{Y}}_{\mathrm{si}}/{\mathrm{Y}}_{\mathrm{pi}}\right)\right]/\left[1-\left({\mathrm{Y}}_{\mathrm{ms}}/{\mathrm{Y}}_{\mathrm{mp}}\right)\right] $$

$$ \mathrm{STI}=\left[\left({\mathrm{Y}}_{\mathrm{pi}}\times {\mathrm{Y}}_{\mathrm{si}}\right)/{\left({\mathrm{Y}}_{\mathrm{mp}}\right)}^2\right] $$

where Y_si designates the yield of the ith genotype grown under deficit irrigation, Y_pi designates that of the ith genotype grown under normal irrigation, Y_ms is the yield mean over all genotypes grown under deficit irrigation, and Y_mp is the yield mean over all genotypes grown under normal irrigation.

Genotyping

Genomic DNA was extracted from young leaf tissues according to the modified method of Murray and Thompson [59]. The quality and concentration of extracted DNA were determined by electrophoresis in 1% agarose gel. Genotyping using SRAP markers was performed following the method of Li and Quiros [9]. Among the SRAP markers available, 30 primer combinations were screened by polymerase chain reaction (PCR). PCR reactions were conducted in volumes of 10 μL, using a BIO-RAD thermocycler. Each PCR reaction was consisted of 1.5 μL of DNA, 1 μL of forward primer, 1 μL of reverse primer, 5 μL of master mix (Amplicon), and 1.5 μL of distilled water. In SRAP analysis, samples were exposed to the following thermic profile: the first five cycles were run at 94 °C for 1 min (denaturing), 35 °C for 1 min (annealing), and 72 °C for 1 min (extension). Then the annealing temperature was raised to 50 °C for another 35 cycles, followed by another extension step of 10 min at 72 °C. Electrophoresis on 12% non-denatured polyacrylamide gels was used for separation of amplified products; and then the products were stained by AgNO3 solution [60].

Statistical analyses

Phenotypic data analysis

Data were tested for normality by Kolmogorov–Smirnov test and homogeneity of variance was tested by Bartlett test. Analysis of variance was performed for the normal and water-deficit irrigations separately, based on the split-plot in time (year) model, using Proc GLM of SAS release 9.4 [61]. Components of variance were calculated for each irrigation regimes, according to the method of Steel and Torrie [62]. Broad-sense heritability (h²_b) was estimated for both irrigation regimes as described by Nguyen and Sleper [63]:

$$ {\mathrm{h}}^2={\upsigma^2}_{\mathrm{g}}/\left({\upsigma^2}_{\mathrm{g}}+{\upsigma^2}_{\mathrm{g}\mathrm{y}}/\mathrm{y}+{\upsigma^2}_{\mathrm{g}\mathrm{r}}/\mathrm{r}+{\upsigma^2}_{\mathrm{e}}/\mathrm{r}\mathrm{y}\right) $$

In which, σ²_g is the genotype, σ²_gy is the genotype × year, σ²_rg is the genotype × rep, and σ²_e is the residual variance, y is the number of years, and r is the number of replicates. The level of genetic variation was estimated with the calculating of phenotypic coefficient of variation (PCV) and genetic coefficient of variation (GCV) as follows:

$$ \mathrm{PCV}=\left({\upsigma}_{\mathrm{p}}/\upmu \right)\ 100 $$

$$ \mathrm{GCV}=\left({\upsigma}_{\mathrm{g}}/\upmu \right)\ 100 $$

where σ_p is the standard deviation of the phenotypic variance, σ_g is the standard deviation of the genotypic variance, and μ is the phenotypic mean [64].

Molecular data analysis

Polymorphic SRAP markers were scored as binary data with presence (1) or absence (0). For each of the SRAP markers, the following indices were computed. Polymorphism information content (PIC) was determined according to the formula of Roldán-Ruiz et al. [65]:

$$ {\mathrm{PIC}}_{\mathrm{i}}=\Big({2\mathrm{f}}_{\mathrm{i}}\times \left[1-\mathrm{f}\right] $$

where i is the ith primer, f_i is the frequency of the amplified allele, and (1-f_i) is the frequency of the null allele. Resolving power (RP) was estimated as:

$$ \mathrm{RP}=\sum 1-\left[2\times \left(0.5-{\mathrm{f}}_{\mathrm{i}}\right)\right] $$

Marker index (MI) was calculated according to Powell et al. [66]:

$$ {\mathrm{MI}}_{\mathrm{i}}={\mathrm{PIC}}_{\mathrm{i}}\times {\mathrm{N}}_{\mathrm{i}}\times {\mathrm{b}}_{\mathrm{i}} $$

where N_i is the total bands for the ith primer, and b_i is the percentage polymorphic bands of the ith primer.

Population structure and association analysis

Structure analysis and stratification of the studied population into subpopulations with different genetic structures was done based on SRAP marker data in STRUCTURE software version 2.3.4 [25]. This analysis was performed applying an admixture model, a burn-in of 10,000 iterations followed by 100,000 Monte Carlo Markov Chain (MCMC) replicates. The membership of each genotype was run for the range of genetic clusters (K) from K = 2 to K = 10 with five replications for each K. delta k approach described by Evanno et al. [67] was used to determine the optimum number of sub-populations, using STRUCTURE HARVESTER [68].

Association analysis was run by both GLM and MLM [26] to calculate P-values for marker–trait associations, using TASSEL version 4.2.1 [69]. The phenotypic mean of traits (P-matrix) over 2 years was used to identify significant associations under normal and water deficit irrigations, separately. To correct for population structure in GLM and MLM models, a Q-matrix that was derived from structure analysis (at maximum DK), was applied as a covariate. Moreover, a kinship matrix (K-matrix) was calculated based on the results of marker genotype data using TASSEL version 4.2.1 [69] and was used in MLM. A correction for multiple testing was done with the FDR (false discovery rate) method, using the QVALUE R package [70].