Skip to main content
  • Research article
  • Open access
  • Published:

Multi-locus genome-wide association mapping for spike-related traits in bread wheat (Triticum aestivum L.)



Bread wheat (Triticum aestivum L.) is one of the most important cereal food crops for the global population. Spike-layer uniformity (the consistency of the spike distribution in the vertical space)-related traits (SLURTs) are quantitative and have been shown to directly affect yield potential by modifying the plant architecture. Therefore, these parameters are important breeding targets for wheat improvement. The present study is the first genome-wide association study (GWAS) targeting SLURTs in wheat. In this study, a set of 225 diverse spring wheat accessions were used for multi-locus GWAS to evaluate SLURTs, including the number of spikes per plant (NSPP), spike length (SL), number of spikelets per spike (NSPS), grain weight per spike (GWPS), lowest tiller height (LTH), spike-layer thickness (SLT), spike-layer number (SLN) and spike-layer uniformity (SLU).


In total, 136 significant marker trait associations (MTAs) were identified when the analysis was both performed individually and combined for two environments. Twenty-nine MTAs were detected in environment one, 48 MTAs were discovered in environment two and 59 MTAs were detected using combined data from the two environments. Altogether, 15 significant MTAs were found for five traits in one of the two environments, and four significant MTAs were detected for the two traits, LTH and SLU, in both environments i.e. E1, E2 and also in combined data from the two environments. In total, 279 candidate genes (CGs) were identified, including Chaperone DnaJ, ABC transporter-like, AP2/ERF, SWEET sugar transporter, as well as genes that have previously been associated with wheat spike development, seed development and grain yield.


The MTAs detected through multi-locus GWAS will be useful for improving SLURTs and thus yield in wheat production through marker-assisted and genomic selection.


Bread wheat (Triticum aestivum L.) is one of the most widely grown cereal crops in the world and serves as the main energy source for approximately one-third of the global population [1,2,3]. Several grain yield parameters play important roles in improving wheat yield, including spike number per plant (SNPP), spike length (SL), grain weight per spike (GWPS), fertile spikelet number (FSN), spikelet number (SN), sterile spikelet number (SSN) and grain number per spike (GNPS). These parameters are important breeding targets for wheat yield improvements [4,5,6,7,8]. Understanding the physiological, genetic and developmental basis of spike morphology is significant not only for increasing spikelet number but also for exploiting the fruiting or grain setting productivity of spikelets [9]. Moreover, characteristics of the spike-layer uniformity (SLU) are usually determined by the variation of spike heights between inter-plants and between inter-tillers. These two indexes of the population uniformity are used not only in wheat but also in rice [10,11,12]. All these factors are affected by environmental factors, agronomic management and genotypes [10, 13].

The most important step for the genetic improvement of wheat is conventional wheat breeding, which is based on phenotypic selection [14]. Breeders usually increase wheat yield by modifying the spike number/hectare, grain number/ear, or 1000 grain weight [15, 16]. However, due to the high cost of generating these traits over large germplasm and the labour involved, the improvement of wheat yield using conventional breeding of spike traits is difficult. Furthermore, these traits are significantly influenced by several other factors, such as genotypic factors, environmental factors and the gene-environment interaction (G X E) [16]. Previously, three important approaches were suggested for increasing the spike size in wheat: (i) increasing spikelet number per spike [17,18,19,20]; (ii) increasing the number of florets and/or grains and grain size per spikelet [7, 21,22,23]; and (iii) simultaneously increasing the spikelet and floret/grain number and grain size. However, combined improvement of these traits is difficult through conventional breeding, as they are often negatively linked to each other. Marker-assisted recurrent selection (MARS) and marker-assisted backcrossing (MABC) are good options. Hence, discovering the important single nucleotide polymorphisms (SNPs)/quantitative trait loci (QTLs) associated with spike-related traits is an urgent task for wheat breeding programs.

Association mapping (AM) was combined with linkage disequilibrium (LD) to identify the function of linked markers and genes for disease-associated loci in humans [24]. Currently, it is also widely used in plants, including cereals [25, 26]. Likewise, through GWAS, marker-trait associations (MTAs) have been identified for yield and its related traits [27, 28]. However, limited efforts have been made in researching the genetics of spike-related traits for wheat. Additionally, only a few studies involving multi-locus and multi-trait GWA analyses have been conducted in wheat in general [29, 30]. There are few studies focussing on wheat yield and quality traits using multi-locus GWAS; therefore, more research is needed to overcome the limitations of routine single-locus single-trait GWAS. In recent years, several molecular mapping studies of spike-related traits have been conducted that led to the identification of QTLs in different wheat germplasms [16, 31,32,33,34,35,36,37], as well as multiple genes that have been identified in wheat for spike morphology traits [18,19,20, 22, 38, 39]. To the best of our knowledge, only a single QTL mapping study has been conducted in wheat for SLURTs involving a RIL population [12].

The present study aimed to identify MTAs involving spike-related traits using LD-based multi-locus GWA mapping using SNP markers in an association panel comprising of 225 spring wheat accessions. Some novel aspects have been included in this study and will be discussed below. First, we included the congruence of our results with historical context which reveals the genes identified for plant architecture and thus identifies the possible genes that could be targeted for yield improvements using SLURT traits. We have also highlighted the novel genes identified in this study using the multi-locus GWAS method. Additionally, we determined that these traits are not independent of the grain traits via providing the grain yield trait correlations from our unpublished data. Finally, we included genomic selection to demonstrate the moderate to high levels of efficiency in predicting SLURTs traits from the identified molecular markers, highlighting the pre-breeding relevance of the genetic markers identified. In this study, the wheat association mapping panel was tested for two years in spring-sown conditions (2017–18 and 2018–19), and phenotypic data were collected for eight different SLURTs (Table 1; Additional File 1: Fig. S1). Candidate genes (CGs) underlying some of the MTAs were also identified. The MTAs identified in this study will be useful for marker-assisted selection (MAS) in the development of new desirable wheat varieties with an ideal spike-layer distribution.

Table 1 List of eight studied traits and their abbreviations


Phenotypic variation and statistical analysis

The ANOVA test revealed significant differences in the measured traits with the following sources of variation: genotypes (G), environments (E), G X E interaction and genotype within replications. Except for LTH, all traits displayed significant variation between the factors, indicating environment-specific effects on measuring these traits, the F-values of each trait are given in Table 2. Thus, to capture the effects of each environment, GWAS analysis was performed on each environment individually along with combined analysis. Estimates of broad-sense heritability (H2) of these traits were generally moderate to high and ranged from 39% (SLN) to 91% (GWPS), indicating the robustness of the measured traits (Table 2).

Table 2 Analysis of variance (ANOVA) of eight spike related traits tested for two environments

The pairwise correlation coefficient (r) analysis of the eight spike-related traits in individual environments, as well as the combined data from both environments (E1 and E2) revealed strong correlations among several measured traits (Additional file 1 Fig. S2, S3, S4). Correlation analysis was also performed for spike-related traits and grain yield (GYPP) from data made available from our unpublished study. Fifty-one significant correlation combinations were found between different traits including grain yield in E1, E2 and the combined data, including 19 in E1 (11 positive and 8 negative) and 15 in E2 (8 positive and 7 negative) and 17 in the combined data (10 positive and 7 negatives). The SLU was negatively correlated with SLN and SLT in E1, E2 and the combined data. SLU was also negatively correlated with grain yield. The correlation between the SLN and SLT was positive in E1, E2 and the combined data. Positive correlations (P < 0.05) were observed between SLN and GYPP for E1 (Additional file 1; Fig. S2) and between SLT and GYPP for E1 and the combined data (Additional file 1; Fig. S2, S4).

Population structure and linkage disequilibrium (LD) analysis

The PCA results illustrate the population structure of the studied population (Fig. 1). Accessions originating from Afghanistan clustered together and formed clusters along the PC1 (towards the right). Accessions from Mexico were distributed almost everywhere. Furthermore, accessions from India and China were also grouped on the PC1 axis. LD decay distance ranged from 2 cM to 20 cM in different genomic regions (30).

Fig. 1
figure 1

Principal component analysis (PCA) obtained using polymorphic SNP markers showing the distribution of 225 SWRS wheat genotypes along with the two components. PC1 explained 7% while PC2 explained 4% variance

Marker–trait associations (MTAs)

A total of 136 significant MTAs were detected for the eight traits in the individual environments and the combined data with (−log10 [p-value] ≥ 3). In environment one (E1), 29 significant MTAs were detected for eight traits, with five for GWPS, seven for LTH, four for NSPP, two for each NSPS & SL, one for each SLN & SLT and seven for SLU. The 29 MTAs were distributed on 11 different chromosomes. In environment two (E2), 48 significant MTAs were found for the eight traits, with three for GWPS, five for LTH, six for NSPP, seven for NSPS, two for SL, seven for SLN, eight for SLT and 10 for SLU. These 48 MTAs were distributed on 14 different chromosomes. In the combined data analysis, 59 significant MTAs were identified for the eight traits, with five each for GWPS & LTH, seven each for NSPP & NSPS, five for SL, 13 for SLN, three for SLT and 14 for SLU. These MTAs were mapped on 18 different chromosomes (Table 3). Altogether, 15 common significant MTAs were found for five traits in either of the two environments. Four MTAs (3 associated with LTH and 1 associated with SLU) were detected in both environments i.e. E1, E2 and also in combined data from the two environments. Fourteen MTAs were common between E2 and the combined data (E1 and E2). A single MTA was common between E1 and the combined data (E1 and E2). Details of all these MTAs are presented in Table 4.

Table 3 List of significant SNP markers (−log10 [p-value] ≥ 3) associated with different spike related traits using the individual environment and combined data
Table 4 List of common significant MTAs for a particular trait in Environment 1 (E1), Environment (E2), and combined environments (E1 and E2) at P-value ≤0.0001

Furthermore, 6 multi-trait MTAs were found to be associated with traits such as SLN, SLT, SLU and GWPS. The details of these MTAs are presented in Table 5. After the false discovery rate (FDR) multiple corrections, only 5 MTAs were identified. These MTAs were detected in E1, combined data and were associated with three traits only. Details of these MTAs are presented in Table 6. Manhattan plots were used to display the SNPs associated with traits. The solid horizontal line indicates the cut-off p-value of significant SNPs with the trait (Fig. 2 and Additional file 1; Fig. S5). Similarly, the quantile-quantile (Q-Q) plots of the observed and expected distribution of p-values found nearly linear trends, showing appropriate model fitting for the GWAS test, as shown in Fig. 3 and Additional file 1 and Fig. S6.

Table 5 Details of common significant MTAs associated with more than one trait
Table 6 List of FDR qualified MTAs for three traits in environment 1 and combined data
Fig. 2
figure 2

Manhattan plots obtained using FarmCPU for eight SLURTs (a) GWPS, (b) LTH, (c) NSPP, (d) NSPS, (e) SL, (f) SLN, (g) SLT (h) SLU under combined E1 and E2 conditions

Fig. 3
figure 3

QQ-plots to visualize the deviation of observed p values from expected p values (based on null hypothesis) for eight SLURTs (a) GWPS, (b) LTH, (c) NSPP, (d) NSPS, (e) SL, (f) SLN, (g) SLT, (h) SLU under combined E1 and E2 conditions

Comparison of MTA or MTA groups (QTLs) detected in the current study to historical QTL regions

In comparison with previously known QTL intervals, MTAs or MTA groups (QTLs) detected in the current study were detected either within or less than 50 MB of previously known QTL intervals. MTA was considered novel if it was detected at a different or faraway (> 50 MB) position. Six MTAs were found within the interval of reported QTLs for the same traits (Additional File 1: Fig. S7). These MTAs were associated with three traits: GWPS [SNP_7209; 40, and SNP_109; 41], LTH [SNP_2448; 42 and SNP_12067; 43], and NSPP [SNP_598 and SNP_454; 44]. Six MTAs were also detected in the interval of QTLs reported previously for grain yield. These MTAs were associated with three traits: GWPS (SNP_5621 and SNP_493; 43], NSPS [SNP_4528 and SNP_13896; 43] and SLN [SNP_2980 and SNP_1140; 43]. Twenty-seven MTAs were also found near the reported QTL intervals for the same traits or grain yield. Fifty three MTAs were novel.

Genomic prediction

To ascertain whether spike-related traits could be predicted using markers, we performed genomic selection (GS) on these traits. Across the eight measured traits, moderate to high prediction values were observed, ranging from r2 = 0.43 to 0.76. The lowest and highest predictions were observed for the NSPP and LTH traits, respectively (Fig. 4).

Fig. 4
figure 4

Genomic predictions for the studied traits

Identification of candidate genes (CGs)

In total, 279 high-confidence CGs were identified using 48 significant non-redundant SNPs (out of 94). However, no hit was found from the remaining 46 SNPs. To identify functions, CGs were annotated using gene ontology (GO) based on IWGSC RefSeq v1.0. These CGs encode a variety of proteins that include genes involved in accelerating developmental stages such as embryogenesis, lateral root development, vascular differentiation, flower development, cell division, and elongation and differentiation in wheat and related species (Additional file 2; Table S3). Furthermore, based on a survey of published literature on wheat’s life cycle considering spike-related traits, 17 CGs were selected and presented (Table 7) [40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78]. The remaining 262 CGs were not associated with the traits considered in the present study. However, these 262 CGs were associated with several other related traits and thus cannot be ignored (Additional file 2; Table S3).

Table 7 List of important candidate genes and their function

We further explored the expression analysis of these genes from public repositories. The results of the in-silico gene expression analysis of these CGs are shown in Fig. 5, with related details provided in Table 7. Given that the traits we analysed in the present study are directly related to plant development, differential gene expression data could be used to identify candidate genes. Among the highly expressed differentially expressed genes, TraesCS4A02G446700, which encodes sucrose synthase and is associated with SLU, showed the highest expression in the stem. TraesCS6A02G334200, which encodes alpha-amylase and is associated with SLU, showed high expression in the spike. CGs TraesCS4B02G272400 and TraesCS4A02G097900, which encode ABC transporter-like and heat shock proteins, respectively, were associated with SLU and showed the highest expression in leaves. TraesCS3B02G309800, which encodes AP2/ERF, was associated with SLN and showed the highest expression in leaves. TraesCS5A02G372900, which encodes Chaperone DnaJ and is associated with GWPS, showed the highest expression in grain. The expression pattern of TraesCS1B02G336900 could not be checked, as the corresponding information was not available in the expression database. This CG was associated with SLN, SLT and SLU and encoded CLAVATA3/ESR (CLE)-related protein 25/26, which is known to be involved in the regulation of cell proliferation and differentiation of plant shoots, roots, vasculature, and more. (Table 7).

Fig. 5
figure 5

Heat map of in silico gene expression analysis for spike trait candidate genes identified through RNA-seq expression data from wheat expression browser (


A major impetus in increasing grain yield comes from finding novel ways of improving crop yields. Improvement in spike-related traits is an important research area that has a direct consequence on grain yield improvements. As screening large segregating populations for these traits in a breeding program can be labour intensive and cost-ineffective, improvements through marker-assisted selection (MAS) or genomic selection is a good option that needs prior marker information or preliminary studies [79, 80]. Thus, understanding the genetics of yield trait components, such as spike morphology, holds the key to wheat yield improvement. In this regard, the availability of wheat genome sequences and subsequently the availability of numerous wheat SNP platforms play important roles in enabling the effective development of cultivars by using genetic resources and genome-enabled breeding.

In wheat, most genome-wide association studies for agronomic and grain yield traits are based on single locus single trait (SLST) analysis [28]. However, it is possible that SLST analysis cannot detect all MTAs [81]. To overcome this limitation, the advanced multi-locus GWAS approach “FarmCPU” (fixed and random model circulating probability unification) has been developed [81]. Hence, in this study, we used FarmCPU, which holds more statistical power and better computational efficiency than other available methods [81] for GWAS analysis, such as EMMA (efficient mixed-model association) [82], EMMAX (EMMA eXpedited) [83] and GEMMA (genome-wide efficient mixed-model association) [84]. FarmCPU eliminates confounding problems arising due to population structure, kinship and multiple testing correction [81, 85,86,87,88,89]. It uses both a fixed-effect model (FEM), a test marker using pseudoQTNs as covariates and a random effect model (REM); tests estimate pseudoquantitative trait nucleotides (QTNs) iteratively [81, 90]. In the FarmCPU model, population structure (PCA) is a fixed effect, kinship is a random effect [81, 88] and model overfitting is avoided by estimating kinship [89]. Given that population structure is well controlled, we did not detect many MTAs for the studied traits. Additionally, Q-Q plots confirmed the suitability of the multi-locus association model used in this study (Fig. 3 and Additional file 1 Fig. S6).

The PCA results obtained in this study suggest that accessions from Mexico harbour greater diversity when compared to India and China, indicating that wheat germplasms from India and China are less diverse. These results are not surprising when considering the large breeding efforts of CIMMYT to diversify their germplasm.

In the present study, we found a total of 136 MTAs for all eight traits in the individual environment and combined data. These MTAs were distributed over 18 wheat chromosomes. Interestingly, a large number of MTAs were from subgenome A, followed by B, and the least was from subgenome D. Since wheat polyploidization involves historical hybridization of the A, B and D genomes, fewer MTAs from the D-genome indicate a founder effect of the wheat population. In the present study, we observed wide ranges in the values of measured traits, as indicative of their distribution (Additional file 2; Table S1). ANOVA of the phenotypic data for different spike-related traits indicated significant genotypic (G) and G X E interactions for most of the studied traits. More MTAs were found in E2 than E1, which indicates the environmental influence of the studied traits. These environmentally specific MTAs may facilitate breeding for specific environmental conditions [91].

Environmental (meteorological) data from both years exhibited large differences in mean temperature and rainfall, suggesting that uneven distribution of these two parameters (over the crop growing period) can lead to large differences in wheat traits (Additional file 1; Table S2). Hence, variation was significantly different in both environments and supports the consideration of two years as two environments. Interestingly, most promising traits, such as GWPS, NSPS, NSPP, SL, LTH and SLU, displayed high heritability and were consequently less influenced by environmental changes, indicating the robustness of these traits. Nevertheless, we measured traits that display strong environmental influence as spike-related traits could be influenced by many environmental variables that make routine wheat breeding difficult. Similar to the present study, genotype-by-environment (G X E) interactions have also been reported previously for traits such as spikelet per spike, spike length, the number of spikes per plant, etc. [12, 92,93,94]. Zhao et al. [12] also observed high heritability (more than 50%) for NSPP, SL, NSPS, GWPS, LTH and SLU, whereas two traits, SLT and SLN, showed heritability of less than 50% (Table 2). Present and previous research suggests the complexity of the SLURTs as they often interact with the environment, demonstrating the difficulty in genetic improvement of SLU by direct selection through conventional wheat breeding programs. Therefore, it is important to perform a genetic analysis of SLURTs in wheat using molecular breeding programs.

In wheat, previous studies [95, 96] have found that plants with bumpy spike layers have a high yield potential. Therefore, uniform spikes (even spike layers) distributed in limited packed horizontal space may reduce photosynthetic rates and yield in wheat. In contrast, uneven spike-layer distribution exposes large leaf surface areas to sunlight causing higher photosynthesis which in turn results in higher wheat yield. Thus, plants with suitable uneven spike layers may be an ideotype for higher yield in wheat [12]. Negative correlations between yield and SLU observed in the present study further support this view. A negative correlation between yield and SLU was also observed in a biparental mapping study [12]. Similarly, Yao [95] and Hu [96] showed that SLU had a significantly negative correlation with yield potential in wheat, and plants with a slight difference in spike heights among tillers (uneven spike-layer to a certain extent) resulted in higher yield potentials. Grain per spike is a primary component of wheat yield, and spikelet number is a primary factor affecting the grains per spike in wheat. The photoperiod and temperature are the two major environmental factors that control spikelet and floret primordia initiation in wheat [17, 97,98,99,100]. Similarly, different developmental stages of spikes affect spikelet number and grain yield [3, 101, 102]. Additionally, spike length variation also affects grain per spike and thus plays an important role in improving wheat yield [14]. We found an overlap of six MTAs for three traits with earlier studies. Six MTAs for three traits were also found within the reported QTL interval for grain yield (Additional file 1 Fig. S7). Moreover, in the present study, four stable MTAs were found (Table 4). Three MTAs were for LTH and one MTA was for SLU. Out of three MTAs for LTH, SNP_2448 was mapped within the reported QTL interval for tiller number [103] on chromosome 1A (Additional file 1 Fig. S7). In addition, six multi-trait MTAs were also detected (Table 5). These findings will provide useful information for wheat breeding programs.

Several significant correlation combinations were observed in this study, such as LTH positively correlated with SLN and SLT in E1 and the combined; similarly, LTH was negatively correlated with SLU in E1 and the combined analysis. Likewise, SLU was negatively correlated with SLN and SLT in both the environments and in the combined analysis. Conversely, a positive correlation between LTH and SLU was found in rice (Oryza sativa) [11], indicating that the genetic mechanisms controlling SLU in wheat and rice are likely under the control of diverse mechanisms. Furthermore, SLN was positively correlated with SLT in all the data sets (Additional file 1 Fig. S2, S3 & S4), a similar pattern of correlation combination was reported by [12] in wheat.

As illustrated in Fig. S1 (left panel) (Additional File 1: Fig. S1), SLN and SLU were identical (value =1) when the genotype within a plot exhibited uniform heights. However, when the spikes are of different heights (right panel of Fig. S1) the SLN and SLU values are not identical. Thus, the value of SLU is 0.5, but the value of SLN is 2. Hence, SL = SLT, SLU = 1, SLN = 1 when all spikes are of identical heights (identical SL). Nevertheless, in this study, we did not detect identical trait values for these traits, a negative or positive correlation occurred between the three derived traits. As both SLU and SLN are inverse of each other, strong negative correlations, − 0.85***, − 0.84*** and − 0.87***, between SLU and SLN were observed for E1, E2 and the combined, respectively (Additional File 1: Fig. S2, S3, S4).

Thus, spike height (spike length or SL) plays an important role in the SLT calculation. Therefore, if all spikes are of the same height, then SLT = SL. However, if spike height (SL) increases or decreases, SLT also correspondingly increases or decreases, respectively. Positive correlations were observed between SL and SLT, 0.24***, 0.21*** and 0.23***, for E1, E2 and combined, respectively. (Additional File 1: Figs. S2, S3, S4). Similarly, SLT increases or decreases depending on the increase or decrease in SLN (as SLN=SLT/SL) and strong positive correlations were observed between SLN and SLT, 0.94***, 0.88*** and 0.94***, for E1, E2 and combined, respectively (Additional File 1: Fig. S2, S3, S4). Likewise, SLU is inversely related to SLT (SLU=SL/SLT); hence, strong negative correlations, − 0.81***, − 0.76*** and − 0.83***, between SLU and SLT were observed for E1, E2 and the combination of E1 and E2, respectively (Additional File 1: Figs. S2, S3, S4). Hence, calculated traits such as SLN and SLU display strong concordance as several significant MTAs overlapped (Table 5).

In this study, after applying FDR correction, only five MTAs were retained (Table 6). It is well known that multiple test correction steps (such as FDR) in addition to reducing the number of undesirable false-positive associations can lead to the loss or disappearance of true associations (false negatives) [104]. Hence, MTAs that disappear after the application of FDR correction should not be ignored entirely. Mainly when MTAs reported within or near previously reported QTL intervals.

It is difficult to know coinciding hits (with previous studies) mainly due to the usage of different mapping populations, population structure, and LD decay. In this study, a conservative threshold of 50 Mbp was used to locate overlapping QTLs. More research is needed to better understand these genomic regions by identifying and characterizing the candidate genes. The present study detected MTAs within or near previously reported QTL intervals, ensuring the robustness of the analysis. It also suggests the possible usage of these potential MTAs in molecular breeding. Large numbers of MTAs identified in the present study were also subjected to further scrutiny to identify the most important MTAs which could be recommended for marker-assisted recurrent selection (MARS) or marker-assisted selection (MAS). MTAs that fulfilled at least one or more of the following criteria were selected for this purpose: (i) lowest P-value, (ii) qualified FDR multiple correction, (iii) common among two or more traits, (iv) stable (identified in both the environment and combined data) and (v) detected in earlier studies (including both interval mapping and GWAS). A total of 16 MTAs were found using the above criteria for different traits (Table 8). One such MTA is SNP_10924. This MTA was associated with SLU, present in both E1 and E2 and the combined environments and surpassed the FDR threshold. Furthermore, the MTAs observed in this study could be grouped into QTL regions considering chromosome-wise linkage disequilibrium (LD). Based on this, five groups were identified involving 14 MTAs (12 SNPs) for three different environments based on LD decay distance on particular chromosomes (Table 9). The distance among the remaining MTAs was greater than the LD decay distance for a particular chromosome. The following MTA groups were observed. (i) In E1, a single group was observed of two SNPs (SNP_1870 and SNP_10789) mapped on the 5A chromosome and associated with SLU. (ii) In E2, two groups were observed. The first group consisted of two SNPs (SNP_4799 and SNP_3408) on chromosome 5B and was associated with SLN and SLU, the second group consisted of four SNPs (SNP_5327, SNP_10576, SNP_8090 and SNP_5909) mapped on the 6D chromosome and associated with SLU. (iii) In the combined data (E1 and E2), two groups were also observed, the first group consisted of two SNPs (SNP_1852 and SNP_2933) mapped on the 6A chromosome and associated with LTH and the second group consisted of two SNPs (SNP_4799 and SNP_3408) mapped on chromosome 5B and associated with SLT. Two SNPs (SNP_4799 and SNP_3408) were observed both for E2 and the combined data.

Table 8 Summary of most important MTAs for MARS or MAS
Table 9 Details of MTAs groups (QTLs) observed on the basis of chromosome-wise LD decay distance

Given that the traits described in this study are environmentally sensitive, we still obtained high genomic prediction values of most of the traits based on the markers. This suggests the merits of genomic selection in breeding for SLURT traits and could be implemented for breeding when labour-intensive phenotyping is not feasible.

An effort was also made to identify CGs using high confidence MTAs, resulting in the identification of a total of 279 CGs, suggesting their possible involvement in different biological processes (Additional file 2; Table S3). Important putative CGs that played important roles in different biological pathways in wheat were selected (Table 7). Domains such as MADS boxes are known to play important roles in reproductive development, regulation, inflorescence architecture, flowering time control, floral organ identity determination and seed development [69, 70]. Likewise, the transcription factor GRAS was also detected, which is well known for its involvement in basic metabolic processes such as photosynthesis, plant growth, senescence and more [72, 73].


The present study is the first report on MTAs for SLURTs in multiple environments using multi-locus GWAS. This study provides insights into the genomic regions of SLURTs. The MTAs identified herein may facilitate breeding new wheat varieties with logically applied spike-layer distribution through MAS. The identified MTAs will be useful for MAS after further validation through suitable mapping populations segregating for these significant MTAs. For high-throughput SNP genotyping, a suitable KASP (Kompetitive Allele-Specific PCR) assay may also be developed using SNP tags of desirable alleles identified through GWAS. CGs identified in the present study should be further validated through functional genomics approaches and may be utilized for developing CG-based markers for their further utilization in CG-based association mapping. In conclusion, we suggest that early-stage selection during breeding programs may focus on the aspect of the suitable vertical spatial distribution of spikes to target higher yield potential.


Wheat association mapping panel and SNP genotyping

The association mapping panel used consisted of 225 diverse wheat genotypes (Additional file 1; Table S1), which is a subset of the original spring wheat reference set (SWRS) of 330 wheat genotypes [109]. Seed material was obtained from the International Maize and Wheat Improvement Centre (CIMMYT), Mexico [109]. Genotyping was performed by CIMMYT through outsourcing using diversity array technology (DArT) in combination with next-generation sequencing platforms known as DArT-seq [109]. Out of the 17,937 SNP markers made available for the original set of 330 genotypes, 9627 (53.67%) were mapped on the genetic map. A total of 5582 markers were finally retained after further filtering the 9627 genetically mapped markers for a minor allele frequency of 5% and missingness (30%).

Field experiment

The association panel of 225 diverse spring wheat genotypes was raised in a simple lattice design with two replications for two years (2017–18 and 2018–19) at the Agriculture Research Farm, Ch. Charan Singh University, Meerut, UP, India; 28.984644°N and 77.705956°E.

Lattice designs are very useful incomplete block designs for plant breeders, as these designs have more flexibility in choosing the number of replications depending upon the availability of resources [110]. The lattice design can be a simple lattice if the design has two replications of the treatments [110]. The replication number and the treatments are flexible in a simple lattice design and thus useful for testing a large number of treatments [111]. Several recent GWAS in crops, including wheat and rice, have used the alpha lattice design with two replications [112, 113].

In this study, each year was treated as one environment, thus making two environments (E1, E2). Data from both years were significantly different for most of the environmental parameters after applying a t-test on both year’s meteorological data (Additional file 1; Table S2). Each genotype was represented by a plot of 3 rows of 1.5 m each, with a row-to-row distance of 0.25 m. The total number of blocks were 15, with each block containing 45 rows, i.e., three rows of each genotype. Five plants per row were used for data recording. First, second, third, fourth and fifth irrigations were conducted 21, 45, 60, 90 and 100 days after sowing. Standard field management practices (i.e., 200 kg/ha fertilizer; N:P:K = 8:8:8) were followed for both environments.


The data on each of the 225 genotypes were recorded for the spike traits according to the following procedure: Number of spikes per plant (NSPP); total number of normal spikes were counted from five randomly selected plants of each plot, out of which five randomly selected spikes were used to calculate other different parameters such as (a) the number of spikelets per spike (NSPS) (spikelet number was counted from the basal sterile spikelet to the top fertile spikelet), (b) spike length (SL) (SL was measured from the base of the spike to the tip, excluding awns), (c) grain weight per spike (GWPS) (after threshing and cleaning of spikes, the grain weight was calculated in grams). Finally, the mean value over the five spikes of each genotype was considered as a final trait value. The lowest tiller height in cm (LTH) was measured from ground level to the tip of the spike (excluding awns), plant height (PH) was measured from ground level to the tip of the spike (excluding awns), and LTH and PH were calculated from a plant selected for NSPP. Similarly, spike-layer thickness (SLT) was calculated from the three traits viz.; PH, LTH and SL (SLT = PH–LTH + SL), spike-layer number (SLN) was calculated as SLN=SLT/SL, and spike-layer uniformity (SLU) was calculated as SLU=SL/SLT. The biological meaning of these traits is presented in Table 1. Measurements for SLURTs were done following [12].

To ascertain the relationship between the spike traits and grain yield per plot (GYPP; in grams), average grain yield data was also measured (size of plots) for two years (2017–18 and 2018–19). The total grain yield of a per m2 plot was weighed in grams. Since the data were considered in another unpublished study on the yield traits, only pairwise Pearson correlation analysis was presented.

Statistical analysis

Descriptive statistics of the measured phenotypes, including the mean, range, standard error and coefficients of variation (CV), were estimated using SPSS Inc. 2008 [114]. Replicated data were used to calculate ANOVA and heritability in a particular year/environment. The mean data of both replications of the individual year/environment of each trait were used to conduct GWAS and correlation analyses. (Additional file 2; Table S2). For the combined analysis over the two environments, least-square means were calculated in the R package “Emmeans” version 1.4.8 ( Pairwise correlation analysis (Pearson’s method) was performed in the R package Performance Analytics among eight traits as well as grain yield. Likewise, the Agricolae R package was used for the estimation of ANOVA using the additive main effects and multiplicative interactions (AMMI) model. Broad sense heritability (H2) was calculated according to [115] from the genotypic variance (σ2g) and phenotypic variance (σ2p) using Microsoft Excel 2010. F-values were calculated for each trait using the environment and replication varieties that were used for the combined analysis of the traits.

Population structure and linkage disequilibrium (LD) analyses

To understand the population structure, PCA was performed from a set of polymorphic SNPs implemented in TASSEL version 5.0 and displayed using [116]. To compute pairwise LD between the markers, TASSEL version 5.2.54 was used. Chromosome-wise, as well as whole-genome LD decay, was conducted in TASSEL 5 from a set of genetically anchored markers, and decay distance (in cM) was plotted. The LD decay of the panel was reported in detail in our previous publication [30].

Marker trait association (MTA) analysis

GWAS was performed using 5582 genetically mapped markers through FarmCPU implemented in GAPIT version 3 [81, 117]. Manhattan and QQ plots were generated using R packages, qqman version 0.1.4 [118]. The P-value threshold was set as 0.001 (−log10 [p-value] =3.0). QQ plots were used to examine model fitting which accounts for population structure. False discovery rate (FDR) correction criteria were also applied to all the identified MTAs with corrected p-values < 0.05 to remove false-positive associations.

Comparison of MTAs detected to historical QTL regions

The MTAs identified were also compared with previously reported QTLs/MTAs for the same traits and grain yield. The physical positions of all MTAs detected in the present study and previously reported QTLs/MTAs were identified through Ensembl Plants [version 50;]. The results were compared based on the physical position of the particular marker on a particular chromosome for the same trait or grain yield. Representative chromosomal maps were prepared for the comparative presentation of all MTAs detected in the present study and previously reported QTLs/MTAs for each trait on an individual chromosome using Map chart software [119].

Genomic selection (GS)

The GS method cBLUP implemented in GAPIT version 3 was used to obtain GS values across the measured traits. Graphics of GS were generated in the R statistical environment using the base package and functions (R Core Team,

Putative candidate gene (CG) identification

A total of 136 significant MTAs were identified in E1, E2 and the combined data, out of which 94 unique MTAs were subjected to the identification of CGs for all the study traits (Additional file 2; Table S3). Candidate genes for these MTAs were identified by aligning the related GBS sequences to wheat genome assembly IWGSC1.0 (IWGSC, 2018), which is hosted on the Ensembl database ( High-confidence annotated genes were retrieved from a chromosome-wise LD decay distance window for each identified MTA. This window size varies for each MTA depending on the local LD decay for each chromosome. The gene ontology (GO) annotation information of these candidate genes (CGs) was extracted from the IWGSC website ( In silico gene expression analysis was also conducted to identify CGs using RNA Seq expression data from Wheat Expression Browser ( A heatmap was generated for the presentation of expression data of genes in different stages and tissues considering the selected domain relevant to the studied traits.

Availability of data and materials

The phenotypic datasets created and analysed during the current investigation for the GWAS are accessible in supplemental table (Additional file 2; Table S2). The genotypic dataset used in this work is accessible from the corresponding author upon reasonable request for non-commercial purposes. As the dataset is subset of the larger panel of accessions from CIMMYT that are genotyped, the corresponding genotypic information can also be obtained from the Nature communication paper [109]; DOI:



Spike-layer uniformity related trait


Genome-wide association study


Association mapping


linkage disequilibrium


Number of spikes per plant


Spike length


Number of spikelets per spike


Grain weight per spike


Lowest tiller height


Spike-layer thickness


Spike- layer number


Spike-layer uniformity


Plant height


Grain yield per plot


Spike number per plant


Fertile spikelet number


Spikelet number


Sterile spikelet number


Grain number per spike


Marker trait association


Single locus single trait


Candidate gene


Single nucleotide polymorphism


Genomic selection


Quantitative trait loci


Diversity arrays technology


Recombinant inbred line


Marker-assisted selection


False discovery rate


Fixed and random model Circulating Probability Unification.


Marker-assisted recurrent selection


Marker-assisted backcrossing


  1. Shiferaw B, Smale M, Braun HJ, Duveiller E, Reynolds M, Muricho G. Crops that feed the world 10. Past successes and future challenges to the role played by wheat in global food security. Food Secur. 2013;5(3):291–317.

    Article  Google Scholar 

  2. Guo J, Shi W, Zhang Z, Cheng J, Sun D, Yu J, et al. Association of yield-related traits in founder genotypes and derivatives of common wheat (Triticum aestivum L.). BMC Plant Biol. 2018;18:1–15.

    Article  CAS  Google Scholar 

  3. Gupta PK, Balyan HS, Sharma S, Kumar R. Genetics of yield, abiotic stress tolerance and biofortification in wheat (Triticum aestivum L.). Theor Appl Genet. 2020;133(5):1569–602.

    Article  PubMed  Google Scholar 

  4. Ijaz U. Smiullah, and Kashif M. genetic study of quantitative traits in spring wheat through generation means analysis. Am J Agric Environ Sci. 2013;13:191–7.

    Google Scholar 

  5. Deng Z, Cui Y, Han Q, Fang W, Li J, Tian J. Discovery of consistent QTLs of wheat spike-related traits under nitrogen treatment at different development stages. Front Plant Sci. 2017;8:1–16.

    Article  Google Scholar 

  6. Guo H, Yan Z, Li X, Xie Y, Xiong H, Liu Y, et al. Development of a high-efficient mutation resource with phenotypic variation in hexaploid winter wheat and identification of novel alleles in the TaAGP.L-B1 gene. Front Plant Sci. 2017;8:1–9.

    Article  Google Scholar 

  7. Liu K, Xu H, Liu G, Guan P, Zhou X, Peng H, et al. QTL mapping of flag leaf-related traits in wheat (Triticum aestivum L.). Theor Appl Genet. 2018;131:839–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Gupta PK, Balyan HS, Sharma S, Kumar R. Biofortification and bioavailability of Zn, Fe and se in wheat: present status and future prospects. Theor Appl Genet. 2021;134(1):1–35.

    Article  CAS  PubMed  Google Scholar 

  9. Slafer GA, Elia M, Savin R, García GA, Terrile II, Ferrante A, et al. Fruiting efficiency: An alternative trait to further rise wheat yield. Food Energy Secur. 2015;4(2):92–109.

    Article  Google Scholar 

  10. Xu ZJ, Huang RD, Li HJ. Difference and correlation of uniformity in rice population among varieties. J Shenyang Agric Univ. 2006;37:137–40.

    Google Scholar 

  11. Ma L, Bao J, Guo L, Zeng D, Li X, Ji Z, et al. Quantitative trait loci for panicle layer uniformity identified in doubled haploid lines of rice in two environments. J Integr Plant Biol. 2009;51(9):818–24.

    Article  CAS  PubMed  Google Scholar 

  12. Zhao C, Zhang N, Wu Y, Sun H, Liu C, Fan X, et al. QTL for spike-layer uniformity and their influence on yield-related traits in wheat. BMC Genet. 2019;20:1–11.

    Article  Google Scholar 

  13. Li JW, Yan QQ. Studies on the relationship between the evenness degree and yield characters of hybrid early rice. Hunan Agric Sci. 2005;1:14–7.

    Google Scholar 

  14. Gao F, Ma D, Yin G, Rasheed A, Dong Y, Xiao Y, et al. Genetic progress in grain yield and physiological traits in Chinese wheat cultivars of southern yellow and Huai Valley since 1950. Crop Sci. 2017;57(2):760–73.

    Article  CAS  Google Scholar 

  15. Ma Z, Zhao D, Zhang C, Zhang Z, Xue S, Lin F, et al. Molecular genetic analysis of five spike-related traits in wheat using RIL and immortalized F2 populations. Mol Gen Genomics. 2007;277(1):31–42.

    Article  CAS  Google Scholar 

  16. Cui F, Ding A, Li J, Zhao C, Wang L, Wang X, et al. QTL detection of seven spike-related traits and their genetic correlations in wheat using two related RIL populations. Euphytica. 2012;186(1):177–92.

    Article  Google Scholar 

  17. Boden SA, Cavanagh C, Cullis BR, Ramm K, Greenwood J, Jean Finnegan E, et al. Ppd-1 is a key regulator of inflorescence architecture and paired spikelet development in wheat. Nat Plants. 2015;1:1–6.

    Article  CAS  Google Scholar 

  18. Dobrovolskaya O, Pont C, Sibout R, Martinek P, Badaeva E, Murat F, et al. Frizzy panicle drives supernumerary spikelets in bread wheat. Plant Physiol. 2015;167(1):189–99.

    Article  CAS  PubMed  Google Scholar 

  19. Poursarebani N, Seidensticker T, Koppolu R, Trautewig C, Gawroński P, Bini F, et al. The genetic basis of composite spike form in barley and ‘miracle-wheat.’ Genetics. 2015;201:155–165.

  20. Dixon LE, Greenwood JR, Bencivenga S, Zhang P, Cockram J, Mellers G, et al. TEOSINTE BRANCHED1 regulates inflorescence architecture and development in bread wheat (Triticum aestivum L.). Plant Cell. 2018;30(3):563–81.

  21. Debernardi JM, Lin H, Chuck G, Faris JD, Dubcovsky J. microRNA172 plays a crucial role in wheat spike morphogenesis and grain threshability. Dev. 2017;144:1966–75.

    CAS  Google Scholar 

  22. Greenwood JR. Wheat inflorescence architecture; 2017.

    Google Scholar 

  23. Sakuma S, Golan G, Guo Z, Ogawa T, Tagiri A, Sugimoto K, et al. Unleashing floret fertility in wheat through the mutation of a homeobox gene. Proc Natl Acad Sci U S A. 2019;116(11):5182–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Bodmer WF, Bailey CJ, Bodmer J, Bussey HJR, Ellis A, Gorman P, et al. Localization of the gene for familial adenomatous polyposis on chromosome 5. Nature. 1988;328:614–6.

    Article  Google Scholar 

  25. Zhu C, Gore M, Buckler ES, Yu J. Status and prospects of association mapping in plants. Plant Genome. 2008;1:5–20.

    Article  CAS  Google Scholar 

  26. Alqudah AM, Sallam A, Baenziger PS, Börner A. GWAS: fast-forwarding gene identification and characterization in temperate cereals: lessons from barley – a review. J Adv Res. 2020;22:119–35.

    Article  PubMed  Google Scholar 

  27. MacCaferri M, Sanguineti MC, Demontis A, El-Ahmed A, Garcia Del Moral L, Maalouf F, et al. Association mapping in durum wheat grown across a broad range of water regimes. J Exp Bot. 2011;62(2):409–38.

    Article  CAS  PubMed  Google Scholar 

  28. Garcia M, Eckermann P, Haefele S, Satija S, Sznajder B, Timmins A, et al. Genome-wide association mapping of grain yield in a diverse collection of spring wheat (Triticum aestivum L.) evaluated in southern Australia. PLoS One. 2019;14:1–19.

    Article  CAS  Google Scholar 

  29. Jaiswal V, Gahlaut V, Meher PK, Mir RR, Jaiswal JP, Rao AR, et al. Genome wide single locus single trait, multi-locus and multi-trait association mapping for some important agronomic traits in common wheat (T. aestivum L.). PLoS One. 2016;11:1–25.

    Article  CAS  Google Scholar 

  30. Kumar J, Saripalli G, Gahlaut V, Goel N, Meher PK, Mishra KK, et al. Genetics of Fe, Zn, β-carotene, GPC and yield traits in bread wheat (Triticum aestivum L.) using multi-locus and multi-traits GWAS. Euphytica. 2018;214:1–17.

    Article  Google Scholar 

  31. Yuan A, Cai C, Lou X, Gao M. Analysis on the genetic model of spike length of wheat main axic. J Luoyang Agric Coll. 1997;17:19–22.

    Google Scholar 

  32. Deng XF, Zhou YH, Yang RW, Ding CB, Zhang L, Zhang HQ, et al. Chromosomal location of genes for spike length in dwarfing polish wheat by monosomic analysis. J Sichuan Agric Univ. 2005;23:12–4.

    Google Scholar 

  33. Marza F, Bai GH, Carver BF, Zhou WC. Quantitative trait loci for yield and related traits in the wheat population Ning7840 x Clark. Theor Appl Genet. 2006;112(4):688–98.

    Article  CAS  PubMed  Google Scholar 

  34. Kumar N, Kulwal PL, Balyan HS, Gupta PK. QTL mapping for yield and yield contributing traits in two mapping populations of bread wheat. Mol Breed. 2007;19(2):163–77.

    Article  Google Scholar 

  35. Li S, Jia J, Wei X, Zhang X, Li L, Chen H, et al. A intervarietal genetic map and QTL analysis for yield traits in wheat. Mol Breed. 2007;20(2):167–78.

    Article  CAS  Google Scholar 

  36. Lu X, Zhang JP, Wang HJ, Yang XM, Li XQ, Li LH. Genetic analysis and QTL mapping of wheat spike traits in a derivative line 3558-2 from wheat x Agronpyron cristatum offspring. J plant Genet Resour. 2011;12:86–91.

    CAS  Google Scholar 

  37. Wang J, Liu W, Wang H, Li L, Wu J, Yang X, et al. QTL mapping of yield-related traits in the wheat germplasm 3228. Euphytica. 2011;177(2):277–92.

    Article  Google Scholar 

  38. Simons KJ, Fellers JP, Trick HN, Zhang Z, Tai YS, Gill BS, et al. Molecular characterization of the major wheat domestication gene Q. Genetics. 2006;172(1):547–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Avni R, Nave M, Barad O, Baruch K, Twardziok SO, Gundlach H, et al. Wild emmer genome architecture and diversity elucidate wheat evolution and domestication. Science. 2017;357(6346):93–7.

    Article  CAS  PubMed  Google Scholar 

  40. Kumar A, Sharma S, Chunduri V, Kaur A, Kaur S, Malhotra N, et al. Genome-wide identification and characterization of heat shock protein family reveals role in development and stress conditions in Triticum aestivum L. Sci Rep. 2020;10:1–12.

    CAS  Google Scholar 

  41. Verslues PE, Zhu JK. Before and beyond ABA: upstream sensing and internal signals that determine ABA accumulation and response under abiotic stress. Biochem Soc Trans. 2005;33(2):375–9.

    Article  CAS  PubMed  Google Scholar 

  42. Finkelstein RR. Studies of abscisic acid perception finally flower. Plant Cell. 2006;18(4):786–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Fujii H, Verslues PE, Zhu JK. Identification of two protein kinases required for abscisic acid regulation of seed germination, root growth, and gene expression in Arabidopsis. Plant Cell. 2007;19(2):485–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. An J, Li Q, Yang J, Zhang G, Zhao Z, Wu Y, et al. Wheat F-box protein TaFBA1 positively regulates plant drought tolerance but negatively regulates stomatal closure. Front Plant Sci. 2019;10:1–20.

    Article  Google Scholar 

  45. Hooykaas PJ, Hall MA, Libbenga KR. Biochemistry and molecular biology of plant hormones: Elsevier; 1999.

  46. Quint M, Gray WM. Auxin signaling. Curr Opin Plant Biol. 2006;9(5):448–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Singla B, Chugh A, Khurana JP, Khurana P. An early auxin-responsive aux/IAA gene from wheat (Triticum aestivum) is induced by epibrassinolide and differentially regulated by light and calcium. J Exp Bot. 2006;57(15):4059–70.

    Article  CAS  PubMed  Google Scholar 

  48. Nicholson P, Chandler E, Draeger RC, Gosman NE, Simpson DR, Thomsett M, et al. Molecular tools to study epidemiology and toxicology of fusarium head blight of cereals. Eur J Plant Pathol. 2003;109(7):691–703.

    Article  CAS  Google Scholar 

  49. Osborne LE, Stein JM. Epidemiology of Fusarium head blight on small-grain cereals. Int J Food Microbiol. 2007;119(1-2):103–8.

    Article  PubMed  Google Scholar 

  50. Gunupuru LR, Arunachalam C, Malla KB, Kahla A, Perochon A, Jia J, et al. A wheat cytochrome P450 enhances both resistance to deoxynivalenol and grain yield. PLoS One. 2018;13:1–17.

    Article  CAS  Google Scholar 

  51. Li Y, Wei K. Comparative functional genomics analysis of cytochrome P450 gene superfamily in wheat and maize. BMC Plant Biol. 2020;20:1–22.

    CAS  Google Scholar 

  52. Lemmens M, Scholz U, Berthiller F, Dall’Asta C, Koutnik A, Schuhmacher R, et al. The ability to detoxify the mycotoxin deoxynivalenol colocalizes with a major quantitative trait locus for fusarium head blight resistance in wheat. Mol Plant-Microbe Interact. 2005;18(12):1318–24.

    Article  CAS  PubMed  Google Scholar 

  53. Walter S, Kahla A, Arunachalam C, Perochon A, Khan MR, Scofield SR, et al. A wheat ABC transporter contributes to both grain formation and mycotoxin tolerance. J Exp Bot. 2015;66(9):2583–93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Gahlaut V, Jaiswal V, Singh S, Balyan HS, Gupta PK. Multi-locus genome wide association mapping for yield and its contributing traits in hexaploid wheat under different water regimes. Sci Rep. 2019;9:1–15.

    Article  CAS  Google Scholar 

  55. Zhang L, Liu P, Wu J, Qiao L, Zhao G, Jia J, et al. Identification of a novel ERF gene, TaERF8, associated with plant height and yield in wheat. BMC Plant Biol. 2020;20:1–12.

  56. Baroja-Fernández E, Muñoz FJ, Saikusa T, Rodríguez-López M, Akazawa T, Pozueta-Romero J. Sucrose synthase catalyzes the de novo production of ADPglucose linked to starch biosynthesis in heterotrophic tissues of plants. Plant Cell Physiol. 2003;44(5):500–9.

    Article  PubMed  Google Scholar 

  57. Keeling PL, Myers AM. Biochemistry and genetics of starch synthesis. Annu Rev Food Sci Technol. 2010;1(1):271–303.

    Article  CAS  PubMed  Google Scholar 

  58. Kim W, Johnson JW, Graybosch RA, Gaines CS. Physicochemical properties and end-use quality of wheat starch as a function of waxy protein alleles. J Cereal Sci. 2003;37(2):195–204.

    Article  CAS  Google Scholar 

  59. Volpicella M, Fanizza I, Leoni C, Gadaleta A, Nigro D, Gattulli B, et al. Identification and characterization of the sucrose synthase 2 gene (Sus2) in durum wheat. Front Plant Sci. 2016;7:1–10.

    Article  Google Scholar 

  60. Cohen P. Protein phosphorylation and hormone action. Proc R Soc B Biol Sci. 1988;234(1275):115–44.

    Article  CAS  Google Scholar 

  61. Schroeder RY, Zhu A, Eubel H, Dahncke K, Witte CP. The ribokinases of Arabidopsis thaliana and Saccharomyces cerevisiae are required for ribose recycling from nucleotide catabolism, which in plants is not essential to survive prolonged dark stress. New Phytol. 2018;217(1):233–44.

    Article  CAS  PubMed  Google Scholar 

  62. Roy S, Vega MV, Harmer NJ. Carbohydrate kinases: a conserved mechanism across differing folds. Catalysts. 2019;9:1–19.

    Article  CAS  Google Scholar 

  63. Ma QJ, Sun MH, Lu J, Zhu XP, Gao WS, Hao YJ. Ectopic expression of apple MdSUT2 gene influences development and abiotic stress resistance in tomato. Sci Hortic. 2017;220:259–66.

    Article  CAS  Google Scholar 

  64. Gautam T, Saripalli G, Gahlaut V, Kumar A, Sharma PK, Balyan HS, et al. Further studies on sugar transporter (SWEET) genes in wheat (Triticum aestivum L.). Mol Biol Rep. 2019;46(2):2327–53.

    Article  CAS  PubMed  Google Scholar 

  65. Tan HT, Shirley NJ, Singh RR, Henderson M, Dhugga KS, Mayo GM, et al. Powerful regulatory systems and post-transcriptional gene silencing resist increases in cellulose content in cell walls of barley. BMC Plant Biol. 2015;15:1–16.

    Article  CAS  Google Scholar 

  66. Kaur S, Dhugga KS, Gill K, Singh J. Novel structural and functional motifs in cellulose synthase (CesA) genes of bread wheat (Triticum aestivum, L.). PLoS One. 2016;11:1–18.

    Article  CAS  Google Scholar 

  67. Galinousky D, Padvitski T, Mokshina N, Gorshkov O, Khotyleva L, Gorshkova T, et al. Expression of cellulose synthase-like genes in two phenotypically distinct flax (Linum usitatissimum L.) subspecies. Genet Resour Crop Evol. 2020;67(7):1821–37.

    Article  CAS  Google Scholar 

  68. Zhu JK. Abiotic stress signaling and responses in plants. Cell. 2016;167(2):313–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Richards RA. Selectable traits to increase crop photosynthesis and yield of grain crops. J Exp Bot. 2000;51(suppl_1):447–58.

    Article  CAS  PubMed  Google Scholar 

  70. Schilling S, Pan S, Kennedy A, Melzer R. MADS-box genes and crop domestication: the jack of all traits. J Exp Bot. 2018;69(7):1447–69.

    Article  CAS  PubMed  Google Scholar 

  71. Liu Y, Zhang H, Han J, Jiang S, Geng X, Xue D, et al. Functional assessment of hydrophilic domains of late embryogenesis abundant proteins from distant organisms. Microb Biotechnol. 2019;12(4):752–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Chen K, Li H, Chen Y, Zheng Q, Li B, Li Z. TaSCL14, a novel wheat (Triticum aestivum L.) GRAS gene, regulates plant growth, photosynthesis, tolerance to photooxidative stress, and senescence. J Genet Genomics. 2015;42(1):21–32.

    Article  CAS  PubMed  Google Scholar 

  73. Niu X, Chen S, Li J, Liu Y, Ji W, Li H. Genome-wide identification of GRAS genes in Brachypodium distachyon and functional characterization of BdSLR1 and BdSLRL1. BMC Genomics. 2019;20:1–18.

    Article  Google Scholar 

  74. Li Z, Liu D, Xia Y, Li Z, Niu N, Ma S, et al. Identification and functional analysis of the CLAVATA3/EMBRYO SURROUNDING REGION (CLE) gene family in wheat. Int J Mol Sci. 2019;20:1–11.

    Google Scholar 

  75. Newberry M, Zwart AB, Whan A, Mieog JC, Sun M, Leyne E, et al. Does late maturity alpha-amylase impact wheat baking quality? Front Plant Sci. 2018;9:1–11.

    Article  Google Scholar 

  76. He H, Hoseney RC. Gas retention in bread dough during baking. Cereal Chem. 1991;68:521–5.

    Google Scholar 

  77. Patel MJ, Ng JHY, Hawkins WE, Pitts KF, Chakrabarti-Bell S. Effects of fungal α-amylase on chemically leavened wheat flour doughs. J Cereal Sci. 2012;56(3):644–51.

    Article  CAS  Google Scholar 

  78. Barrera GN, Tadini CC, León AE, Ribotta PD. Use of alpha-amylase and amyloglucosidase combinations to minimize the bread quality problems caused by high levels of damaged starch. J Food Sci Technol. 2016;53(10):3675–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Venske E, dos Santos RS, Busanello C, Gustafson P, Costa de Oliveira A. Bread wheat: a role model for plant domestication and breeding. Hereditas. 2019;156(1):16.

    Article  PubMed  PubMed Central  Google Scholar 

  80. Wang X, Xu Y, Hu Z, Xu C. Genomic selection methods for crop improvement: current status and prospects. Crop J. 2018;6(4):330–40.

    Article  Google Scholar 

  81. Liu X, Huang M, Fan B, Buckler ES, Zhang Z. Iterative usage of fixed and random effect models for powerful and efficient genome-wide association studies. PLoS Genet. 2016;13:1–24.

    Google Scholar 

  82. Kang HM, Zaitlen NA, Wade CM, Kirby A, Heckerman D, Daly MJ, et al. Efficient control of population structure in model organism association mapping. Genetics. 2008;178(3):1709–23.

    Article  PubMed  PubMed Central  Google Scholar 

  83. Kang HM, Sul JH, Service SK, Zaitlen NA, Kong SY, Freimer NB, et al. Variance component model to account for sample structure in genome-wide association studies. Nat Genet. 2010;42(4):348–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for association studies. Nat Genet. 2012;44(7):821–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Kaur S, Zhang X, Mohan A, Dong H, Vikram P, Singh S, et al. Genome-Wide Association Study reveals novel genes associated with culm cellulose content in bread wheat (Triticum aestivum, L.). Front. Plant Sci. 2017;8:1913.

    Article  PubMed  PubMed Central  Google Scholar 

  86. Gyawali A, Shrestha V, Guill KE, Flint-Garcia S, Beissinger TM. Single-plant GWAS coupled with bulk segregant analysis allows rapid identification and corroboration of plant-height candidate SNPs. BMC Plant Biol. 2019;19(1):412.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. Ward BP, Brown-Guedira G, Kolb FL, van Sanford DA, Tyagi P, Sneller CH, et al. Genome-wide association studies for yield-related traits in soft red winter wheat grown in Virginia. PLoS One. 2019;14(2):e0208217.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. Zhang Y, Wan J, He L, Lan H, Li L. Genome-wide association analysis of plant height using the maize F1 population. Plants (Basel). 2019;8:432.

    Article  CAS  Google Scholar 

  89. Muhammad A, Li J, Hu W, Yu J, Khan SU, Khan MHU, et al. Uncovering genomic regions controlling plant architectural traits in hexaploid wheat using different GWAS models. Sci Rep. 2021;11(1):6767.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Li C, Fu Y, Sun R, Wang Y, Wang Q. Single-locus and multi-locus genome-wide association studies in the genetic dissection of fiber quality traits in upland cotton (Gossypium hirsutum L.). Front Plant Sci. 2018;9:1–16.

    Google Scholar 

  91. Bao B, Chao H, Wang H, Zhao W, Zhang L, Raboanatahiry N, et al. Stable, environmental specific and novel QTL identification as well as genetic dissection of fatty acid metabolism in Brassica napus. Front Plant Sci. 2018;9:1018.

    Article  PubMed  PubMed Central  Google Scholar 

  92. Ding AM, Li J, Cui F, Zhao CH, Ma HY, Wang HG. Mapping QTLs for yield related traits using two associated RIL populations of wheat. Acta Agron Sin. 2011;37:1511–24.

    CAS  Google Scholar 

  93. Zhou Y, Conway B, Miller D, Marshall D, Cooper A, Murphy P, et al. Quantitative trait loci mapping for spike characteristics in hexaploid wheat. Plant Genome. 2017;10:1–15.

    Article  CAS  Google Scholar 

  94. Deng M, Wu F, Zhou W, Li J, Shi H, Wang Z, et al. Mapping of QTL for total spikelet number per spike on chromosome 2D in wheat using a high-density genetic map. Genet Mol Biol. 2019;42(3):603–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Yao WC. Studies on inheritance of spike layer uniformity of wheat. Seed. 2000;109:19–21.

    Google Scholar 

  96. Hu YJ. The difference of the spike-layer architecture and its relation to yield in winter wheat cultivars. Seed. 2001;119:19–21.

    Google Scholar 

  97. Rahman MS, Wilson JH. Determination of spikelet number in wheat. I. Effect of varying photoperiod on ear development. Aust J Agric Res. 1977;28(2):265–74.

    Article  Google Scholar 

  98. Rahman MS, Wilson JH. Determination of spikelet number in wheat. III.* effect of varying temperature on ear development. Aust J Agric Res. 1978;29(3):459–67.

    Article  Google Scholar 

  99. Baker CK, Gallagher JN. The development of winter wheat in the field 1. Relation between apical development and plant morphology within and between seasons. J Agric Sci. 1983;101(2):327–35.

    Article  Google Scholar 

  100. Steinfort U, Fukai S, Trevaskis B, Glassop D, Chan A, Dreccer MF. Vernalisation and photoperiod sensitivity in wheat: the response of floret fertility and grain number is affected by vernalisation status. F Crop Res. 2017;203:243–55.

    Article  Google Scholar 

  101. Rawson HM. Spikelet number, its control and relation to yield per ear in wheat. Aust J Biol Sci. 1970;23(1):1–16.

    Article  Google Scholar 

  102. Baker CK, Gallagher JN. The development of winter wheat in the field: 2. The control of primordium initiation rate by temperature and photoperiod. J Agric Sci. 1983;101(2):337–44.

    Article  Google Scholar 

  103. Ren T, Hu Y, Tang Y, Li C, Yan B, Ren Z, et al. Utilization of a Wheat55K SNP array for mapping of major QTL for temporal expression of the tiller number. Front Plant Sci. 2018;9:1–12.

    Article  Google Scholar 

  104. Kaler AS, Purcell LC. Estimation of a significance threshold for genome-wide association studies. BMC Genomics. 2019;20:1–8.

    Article  CAS  Google Scholar 

  105. Mason RE, Mondal S, Beecher FW, Pacheco A, Jampala B, Ibrahim AMH, et al. QTL associated with heat susceptibility index in wheat (Triticum aestivum L.) under short-term reproductive stage heat stress. Euphytica. 2010;174:423–36.

    Article  Google Scholar 

  106. Li F, Wen W, Liu J, Zhang Y, Cao S, He Z, et al. Genetic architecture of grain yield in bread wheat based on genome-wide association studies. BMC Plant Biol. 2019;19:1–19.

    Google Scholar 

  107. Huang XQ, Cöster H, Ganal MW, Röder MS. Advanced backcross QTL analysis for the identification of quantitative trait loci alleles from wild relatives of wheat (Triticum aestivum L.). Theor Appl Genet. 2003;106:1379–89.

    Article  CAS  PubMed  Google Scholar 

  108. Guan P, Lu L, Jia L, Kabir MR, Zhang J, Lan T, et al. Global QTL analysis identifies genomic regions on chromosomes 4A and 4B harboring stable loci for yield-related traits across different environments in wheat (Triticum aestivum L.). Front Plant Sci. 2018;9:1–18.

    Article  Google Scholar 

  109. Sansaloni C, Franco J, Santos B, Percival-Alwyn L, Singh S, Petroli C, et al. Diversity analysis of 80,000 wheat accessions reveals consequences and opportunities of selection footprints. Nat Commun. 2020;11(1):4572.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  110. Singh P, Bhatia D. Incomplete block designs for plant breeding experiments. Agric Res J. 2017;54(4):607–11.

    Article  Google Scholar 

  111. Patterson HD, Williams ER. A new class of resolvable block designs. Biometrika. 1976;63(1):83–92.

    Article  Google Scholar 

  112. Rahimi Y, Bihamta MR, Taleei A, Alipour H, Ingvarsson PK. Genome-wide association study of agronomic traits in bread wheat reveals novel putative alleles for future breeding programs. BMC Plant Biol. 2019;19:1–19.

    Article  CAS  Google Scholar 

  113. Bhandari A, Sandhu N, Bartholome J, Cao-Hamadoun TV, Ahmadi N, Kumari N, et al. Genome-wide association study for yield and yield related traits under reproductive stage drought in a diverse indica-aus Rice panel. Rice. 2020;13:1–22.

    Article  Google Scholar 

  114. SPSS Inc. Released. 2008. SPSS statistics for windows, version 17.0. Chicago: SPSS Inc.

  115. Allard RW. Principles of plant breeding. 2nd ed: John Wiley & Sons; 1999.

  116. R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R foundation for statistical computing; 2013.

    Google Scholar 

  117. Lipka AE, Tian F, Wang Q, Peiffer J, Li M, Bradbury PJ, et al. GAPIT: genome association and prediction integrated tool. Bioinformatics. 2012;28(18):2397–9.

    Article  CAS  PubMed  Google Scholar 

  118. Turner SD. qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. J Open Source Softw. 2018;3:731.

    Article  Google Scholar 

  119. Voorrips RE. Mapchart: software for the graphical presentation of linkage maps and QTLs. J Hered. 2002;93(1):77–8.

    Article  CAS  PubMed  Google Scholar 

Download references


The authors are thankful to CIMMYT, Mexico for providing the research material. Thanks are due to the Department of Biotechnology (DBT), Govt of India for providing funds in the form of research projects awarded to Shailendra Sharma (SS1). Authors are also thankful to Dr. PK Gupta and Dr. HS Balyan (emeritus professors at CCSU, Meerut) for their constant support. The authors are also thankful to Ch. Charan Singh University, Meerut for providing laboratory and field facilities. The authors declare no conflicts of interest.


This work was supported by the departmental funds received from Ch. Charan Singh University, Meerut to SS1.

Author information

Authors and Affiliations



SS1 (Shailendra Sharma) conceived the research idea and drafted the manuscript. PM, JK, SS2 (Shiveta Sharma) and SS1 performed field experiments. PM, JK, RS and SS2 performed data analysis and also contributed to the writing and revision of the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Shailendra Sharma.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1:

List of 225 spring wheat reference set genotypes. Table S2: Comparative analysis of environmental factors in two different environments/years using t-test. Fig. S1. The biological meaning of spike layer uniformity (the consistency of the spike distribution in the vertical space) related traits (SLURTs) shown with wheat plants having similar tillers height (a) and different tiller height (b). In (a) only one spike layer can be seen from the vertical perspective (SLN = 1) as spikes per plant have consistent vertical distribution (SLU = 1) and the SLT was identical to one spike length (SL) whereas in (b) two spike layer can be seen from the vertical perspective (SLN = 2) as spikes per plant have inconsistent vertical distribution (SLU = 0.5) and the SLT was identical to two spike length. A detailed description of SLURTs is provided in Table 1. Note: In the figure, immature wheat plants are shown just to illustrate the biological meaning of SLURTs. However, data for SLURTs were recorded on mature plants. In the figure scaling is approximate. The figure is based on the publication by Zhao et al. [12]. Fig. S2. Pairwise Pearson correlation among the eight SLURTs in E1 and grain yield (GYPP). A single asterisk (*) represents 0.05 level of significance. A double asterisk (**) represents 0.01 level of significance and a triple asterisk (***) represents 0.001 level of significance. Fig. S3. Pairwise Pearson correlation among the eight SLURTs in E2 and grain yield (GYPP). A single asterisk (*) represents 0.05 level of significance. A double asterisk (**) represents 0.01 level of significance and a triple asterisk (***) represents 0.001 level of significance. Fig. S4. Pairwise Pearson correlation among the eight SLURTs and grain yield (GYPP) for the combined E1 and E2 environments. A single asterisk (*) represents 0.05 level of significance. A double asterisk (**) represents 0.01 level of significance and a triple asterisk (***) represents 0.001 level of significance. Fig. S5. Manhattan plots obtained after using FarmCPU for eight SLURTs (a) GWPS (b) LTH (c) NSPP (d) NSPS (e) SL (f) SLN (g) SLT (h) SLU under E1 (above numbers) and E2 (below numbers) conditions. Numbers correspond to wheat chromosome. Fig. S6. QQ-plots to visualize the deviation of observed p values from expected p values (based on null hypothesis) for eight SLURTs (a) GWPS (b) LTH (c) NSPP (d) NSPS (e) SL (f) SLN (g) SLT (h) SLU under E1 and E2 conditions. Fig. S7 (i-v). Comparison of MTAs or MTA groups, detected in the present study, with historical QTLs/MTAs in different chromosomes; MTAs detected in known flanking regions or locations of QTLs/MTAs for same traits or grain yield are depicted in green color and MTAs found near to flanking regions or locations (less than 50 Mb) of QTLs/MTAs depicted in red colour. Flanking markers are shown with the same colour and with the reference details of the previous study. Novel MTAs detected in the present study are presented without any colour. The corresponding physical distances (Mb) of the QTL/MTA regions on each chromosome were obtained by blasting the flanking sequences of markers, depicted on the right side of each figure, to the Chinese Spring RefSeq v1.0. GY; grain yield, TN; tiller number.

Additional file 2: Table S1.

Descriptive statistics details for all the studied traits. Table S2. Phenotypic data of environment 1 (E1: 2017–18), environment 2 (E2: 2018–19) and mean of E1 and E2. Table S3. List of candidate genes identified using MTAs detected in E1, E2 and combined (E1 and E2).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Malik, P., Kumar, J., Sharma, S. et al. Multi-locus genome-wide association mapping for spike-related traits in bread wheat (Triticum aestivum L.). BMC Genomics 22, 597 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: