Uncovering Genomic Regions Associated With 36 Agro-Morphological Traits in Indian Spring Wheat Using GWAS

Sheoran, Sonia; Jaiswal, Sarika; Kumar, Deepender; Raghav, Nishu; Sharma, Ruchika; Pawar, Sushma; Paul, Surinder; Iquebal, M. A.; Jaiswar, Akanksha; Sharma, Pradeep; Singh, Rajender; Singh, C. P.; Gupta, Arun; Kumar, Neeraj; Angadi, U. B.; Rai, Anil; Singh, G. P.; Kumar, Dinesh; Tiwari, Ratan

doi:10.3389/fpls.2019.00527

ORIGINAL RESEARCH article

Front. Plant Sci., 25 April 2019

Sec. Plant Breeding

Volume 10 - 2019 | https://doi.org/10.3389/fpls.2019.00527

This article is part of the Research Topic Use of Barley and Wheat Reference Sequences: Downstream Applications in Breeding, Gene Isolation, GWAS and Evolution View all 19 articles

Uncovering Genomic Regions Associated With 36 Agro-Morphological Traits in Indian Spring Wheat Using GWAS

$\r\nSonia Sheoran&#x;$ Sonia Sheoran^1†

Sarika Jaiswal^2†

¹ICAR-Indian Institute of Wheat and Barley Research, Karnal, India
²ICAR-Indian Agricultural Statistics Research Institute, New Delhi, India
³Lokbharti-Sanosara Centre, Bhavnagar, India

Wheat genetic improvement by integration of advanced genomic technologies is one way of improving productivity. To facilitate the breeding of economically important traits in wheat, SNP loci and underlying candidate genes associated with the 36 agro-morphological traits were studied in a diverse panel of 404 genotypes. By using Breeders’ 35K Axiom array in a comprehensive genome-wide association study covering 4364.79 cM of the wheat genome and applying a compressed mixed linear model, a total of 146 SNPs (-log₁₀ P ≥ 4) were found associated with 23 traits out of 36 traits studied explaining 3.7–47.0% of phenotypic variance. To reveal this a subset of 260 genotypes was characterized phenotypically for six quantitative traits [days to heading (DTH), days to maturity (DTM), plant height (PH), spike length (SL), awn length (Awn_L), and leaf length (Leaf_L)] under five environments. Gene annotations mined ∼38 putative candidate genes which were confirmed using tissue and stage specific gene expression data from RNA Seq. We observed strong co-localized loci for four traits (glume pubescence, SL, PH, and awn color) on chromosome 1B (24.64 cM) annotated five putative candidate genes. This study led to the discovery of hitherto unreported loci for some less explored traits (such as leaf sheath wax, awn attitude, and glume pubescence) besides the refined chromosomal regions of known loci associated with the traits. This study provides valuable information of the genetic loci and their potential genes underlying the traits such as awn characters which are being considered as important contributors toward yield enhancement.

Introduction

Wheat (Triticum aestivum L.) crop provides one-fifth of total food calories and a quarter of protein in the human diet on daily basis¹. To meet the increasing food demand of growing population, the breeders focused on the varieties having higher yield and yield stability, increased resistance/tolerance to biotic and abiotic stresses. Approximately 10,000 wheat varieties worldwide² including 448 wheat varieties in India (Gupta et al., 2018) have been notified. Agro-morphological characterization of germplasm is fundamental in order to provide information for plant breeding programs. The QTL mapping methods based on bi-parental mapping populations identify the genomic regions with low resolution, whereas, Genome-wide association studies (GWAS), based on linkage disequilibrium (LD), take diverse genetic background into consideration to dissect the genetic architecture of complex traits with high resolution. The GWAS in wheat has started gaining importance in the recent past mainly focusing on yield and yield related traits (Liu et al., 2014; Liu Y. et al., 2017; Sukumaran et al., 2014, 2018; Arruda et al., 2015; Gao et al., 2015; Maccaferri et al., 2015; Arora et al., 2017).

Advances in next generation sequencing technology provided valuable wheat genomic and plant breeding resources including high quality genome data (Brenchley et al., 2012; Jia et al., 2013; International Wheat Genome Sequencing Consortium (IWGSC), 2014; Chapman et al., 2015). Several high throughput SNP arrays viz., 9K (Cavanagh et al., 2013), 90K (Wang, 2014), 820K (Winfield et al., 2016), 660K (Cui et al., 2017), 35K (Allen et al., 2017), and TaBW280K (Rimbert et al., 2018) have been developed and utilized in wheat. These SNP arrays have been successfully used for GWAS in European winter and spring wheat (Zanke et al., 2014), CIMMYT spring wheat (Sukumaran et al., 2014), United States elite wheat breeding genotypes (Lin et al., 2016), a panel of CIMCOG (CIMMYT Mexico core germplasm) Kazakhstan, Russian, and European wheat genotypes (Turuspekov et al., 2017) and Chinese bread wheat cultivars (Sun et al., 2017). A substantial number of novel SNP variants have been identified using 35K SNP Breeders arrays in the Watkins collection of landraces for further improvement of modern elite cultivars (Winfield et al., 2018). From the 820K SNP array (using global selection of germplasm including elite cultivars, landraces, progenitor, and ancestral species of wheat), Breeder’s 35K Axioms array was developed which contains only mapped SNPs that are tailored to be most informative for specific purposes (Wilkinson et al., 2012; Borrill et al., 2015). 35K SNP array holds promise for detecting large scale variation in secondary and tertiary gene pools (Rasheed et al., 2017).

There are several agro-morphological traits which have been studied intensively and for which markers have been identified (Sukumaran et al., 2018). At the same time there are certain less explored traits for instance awn characters which can be considered as an alternative target for the improvement of wheat grain yield through their known functions including photosynthesis and increased water use efficiency (Rebetzke et al., 2016). The present study includes these characters of future importance hitherto not explored much until now. Moreover, most of the agro-morphological traits undertaken in this study are also utilized for the characterization of the genotypes using Distinctiveness Uniformity and Stability (DUS)³.

Some of the studied traits and the associated markers will be of immense importance in future toward developing input use efficient wheat varieties. Amongst them, genomic regions associated with days to heading (DTH) shall enable development of early maturing wheat genotypes to avoid terminal heat stress and allowing an intervening legume crop before rice in the ensuing season (Tewolde et al., 2006; Joshi et al., 2007). For GWAS, size and diversity of the panel plays a significant role as it is suggested that the smaller panel (<384 accessions) and large LD blocks identified in association studies may lead to the identification of false positive associations (Turuspekov et al., 2017). Keeping this in account, a panel of 404 diverse genotypes comprising of indigenous collections, local landraces, released varieties and other improved genotypes was used. The panel also included registered genetic stocks characterized for early maturity, resistance/tolerance to biotic/abiotic stresses, adaptation to different environments, plant architecture etc. (Kundu et al., 2010). These are the major drivers for trait improvement program, using molecular breeding approach. Moreover, India has unique climatic variations rendering wheat germplasm diversity as gold mines. Genotypes adapted to different agro-climatic zones of the country are present in the genotypic panel. Therefore these may be considered as representative of the three major spring wheat growing mega-environments viz., ME 1, ME 4, and ME 5 described by the International Maize and Wheat Improvement Center (CIMMYT) spanning across all the 5 continents⁴. This will allow the breeders to utilize the information in developing genotypes specific to different adaptation conditions.

The present GWA study is the first attempt to undertake large scale evaluation of 35K Axiom array in a diverse panel of 404 Indian wheat genotypes. The 35K Breeders’ Array was selected for present work due to its proven efficacy on panel of 1807 accessions of hexaploid wheat (804 accessions from Watkins Collection, a collection of wheat landraces made by A. E. Watkins in the 1920s and 1930s and 1003 modern and elite accessions) from 32 countries (Winfield et al., 2018). The aim of the study was to (i) identify significant MTAs for 36 agro-morphological traits for future breeding and (ii) mining putative candidate genes underlying the corresponding traits of interest. Furthermore, tissue and growth stage-specific gene expression data was also examined providing support to the detected candidate genes. For complex quantitative traits, the association panel was phenotyped at two locations for 2 years. The outcome of this study could be used to make effective strategies for the development of new varieties coupled with economic traits.

Materials and Methods

Plant Material

A set of 404 bread wheat (T. aestivum L.) genotypes comprising of indigenous collections (91), landraces (8), released varieties (134), genetic stocks (43), and improved genotypes (128) was used to constitute a diverse association panel. The diverse lines were selected on the basis of pedigree, to reduce associations of spurious markers as it provides a buffer against skewness in terms of the environmental effects. Recently, we analyzed trait based diversity analysis using Shannon Index with 16 traits out of 36 traits with a collection of 7,000 diverse germplasm lines (unpublished work). Out of these, 450 genotypes comprising of indigenous collections, landraces, released varieties, improved genotypes and genetic stocks for various traits were selected from 7,000 lines which was downsized to 404 genotypes after eliminating admixtures, duplicates, etc. Comparing the trait based diversity of these selected 404 genotypes using Shannon Index endorsed it as good representative of 7,000 germplasm lines, thereby proving the suitability of 404 genotypes for GWA study (Supplementary Table S1). Seeds of 404 genotypes were obtained from the Germplasm Resource Unit, ICAR-IIWBR (Indian Institute of Wheat and Barley Research), Karnal, Haryana, India, which acts as a nodal center for wheat in the country. Detailed information with pedigree for each genotype is given in Supplementary Table S2.

Field Trials and Phenotyping

The 404 genotypes were evaluated for 30 qualitative characters at the experimental field of ICAR-IIWBR, Karnal during the crop season 2016–2017. A subset of 260 genotypes was phenotyped for six quantitative traits [days to heading (DTH), days to maturity (DTM), spike length (SL), plant height (PH), awn length (Awn_L), and leaf length (Leaf_L)] at three locations viz., Experimental field, Karnal (29°42′N, 77°02′E); Seed Farm, Karnal (29.7138° N, 76.9943°E), Haryana, India and Lokbharti-Sanosara Centre, Bhavnagar (21°46′ N 72°11′ E), Gujarat, India during year 2016–2017. Besides these three environments, an additional environment was taken by phenotyping the subset in crop season 2017–2018 at ICAR-IIWBR, Karnal. Experiment was conducted in two replications following alpha lattice design. To minimize the variations, every genotype was planted with a dibbling tool named IIWBR Dibbler (Sharma et al., 2016) having four rows. The plant to plant distance was 10 cm and row to row distance was maintained at 20 cm. This unique sowing method has helped in avoiding confounding effects of extraneous errors and improved the precision in phenotyping leading to moderate to high estimation of heritability (H²) thereby enhancing the probability of identifying genes of minor effects related to complex traits (Sharma et al., 2016).

Thirty six agro-morphological characters including coleoptiles anthocyanin coloration (C_Col), plant growth habit (PGH), foliage color (Fol_Col), flag leaf anthocyanin coloration of auricle (Aur_Col), flag leaf hairs on auricle (pubescence) (Aur_Pub), flag leaf attitude (Leaf_Att), ear emergence/days to heading (DTH), flag leaf waxiness of sheath (Wax_LS), flag leaf waxiness of blade (Wax_LB), ear waxiness (Wax_Ear), waxiness of peduncle (Wax_Ped), flag leaf length (Leaf_L), flag leaf breath (width) (Leaf_Br), PH, ear shape (Ear_S), ear density (Ear_D), ear (spike) length (SL), awn presence (Awn_P), awns length (Awn_L), awn color (Awn_Col), awn attitude (Awn_Att), outer glume pubescence (Glu_Pub), ear color (Ear_Col), lower glume shoulder width (Sh_Wid), lower glume: shoulder shape (Shl_Sh), beak length (Beak_L), beak shape (Beak_Sh), spike (peduncle) attitude (Ped_Att), grain coloration with phenol (Grn_Ph), grain color (Grn_Col), grain shape (Grn_Sh) grain germ width (Germ_Wid), brush hair length (Brush_L), seed (grain) size (Grn_Size), grain hardness (texture) (Grn_Tex) and DTM were recorded as per guidelines laid out by Protection of Plant Varieties and Farmers’ Right Authority (PPV and FRA, 2011)⁵. Procedure of recording the data for each trait is summarized in Supplementary File S1. Qualitative traits were recorded as binary (presence or absence), ordinal (visual scale of the expression intensity of a characteristic) and nominal (color or shape) (Supplementary Table S1). For association analysis, a total of five environments (E1–E5) were considered namely, E1 – average of two replications at ICAR-IIWBR, Karnal (2016); E2 – average of two replications at Seed Farm, Karnal (2016); E3–average of two locations of Karnal (2016); E4 – average of two replications at ICAR-IIWBR, Karnal (2017) and E5 – average of two replications at Bhavnagar, Gujarat (2016).

Statistical Analysis

Phenotypic data was analyzed using SAS v.9.3 (SAS Institute 2011⁶). Pearson pairwise correlation was calculated for all the traits. Histograms were created in R (R Development Core Team, 2013) using the hist() function. The PROC CORR procedure was employed to calculate correlations among phenotypes. Variance components for the quantitative traits were analyzed using general linear model to detect the effect of genotypes, environment, replication and genotype × environment interaction. All sources of variation were considered as random effects. The broad sense heritability for the traits was estimated by the formula H² = V_G/(V_G + G_E) where V_G and V_E represent estimates of genetic and environmental variance, respectively.

SNP Genotyping and Filtering

Genomic DNA was extracted from 15 days old seedlings according to the CIMMYT Molecular Genetics Manual (Dreisigacker et al., 2013). A Nanodrop 1000 spectrophotometer was used for quantifying DNA at 260 nm absorbance (Biodrop Touch PC+125). The DNA samples were used for genotyping with 35K Axiom^® Wheat Breeder’s Array (Affymetrix UK Ltd., United Kingdom). Quality preprocessing of 35,143 markers obtained from 35K chip was done by using PLINK software (Purcell et al., 2007). Markers with more than 5% missing values, less than 5% minor allele frequency (MAF) and individuals with more than 15% missing SNP calls were removed from the dataset. Markers with no known chromosomal positions, based on high density consensus map generated by using five mapping populations⁷ (Allen et al., 2017), were also removed. Duplicate markers were further removed by R/QTL software (Broman et al., 2003; Arends et al., 2010).

Genetic Diversity and Population Structure Analysis

The basic statistics such as genetic diversity (GD) and polymorphism information content (PIC) was evaluated by using PowerMarker v3.2.5 (Liu and Muse, 2005). The model-based Bayesian cluster analysis program, STRUCTURE v2.3.4 (Hubisz et al., 2009) was used to infer the population structure. A total 100,000 burn-in periods followed by 100,000 Markov Chain Monte Carlo (MCMC) iterations from K = 2 to K = 7 clusters were used to identify the optimal cluster (K). Three independent runs were generated for each K. The results of the analysis were used as input to the Structure Harvester tool (Earl and VonHoldt, 2011) to predict the best K-value based on Evanno method (Evanno et al., 2005). Principal component analysis (PCA) and Neighbor-joining (NJ) tree were created to validate population stratification with the software GAPIT (Lipka et al., 2012) and DARwin v6 (Perrier and Jacquemoud-Collet, 2006), respectively.

Linkage Disequilibrium

For Linkage disequilibrium analysis, r² (squared correlation coefficient) values among markers of all pairs of loci were calculated using PLINK 1.9 tool⁸ (Purcell et al., 2007). Default window size cut off of r² value was used for this analysis. Finally, LD plotting was done for three sub genomes (A, B, and D genomes), on the basis of centiMorgans (cM) distance, using ggplot2 package of R Bioconductor (Wickham et al., 2009). The percentage of marker pairs below the critical LD (r² > 0.02) was also compared in the sub-genomes. Pairwise LD estimates in the region of interest for significantly associated markers were investigated using Haploview 4.2 (Barrett et al., 2005).

Association Analysis

Association analysis was performed using compressed mixed linear model (CMLM) implemented by Genomic Association and Prediction Integrated Tool (GAPIT) in R (Lipka et al., 2012) which took into account a K-PC model (Zhao et al., 2007) where kinship information together with the first three principal components (PC) as covariates were included for GWAS, which further improves statistical power. Kinship matrix was iteratively calculated using the VanRaden method (VanRaden, 2008). The best fit of the model was evaluated on the Q-Q plots generated by the model. A threshold of -log₁₀ P > 4 (-log₁₀P ≥ 4 for quantitative traits) was used to state significant marker trait associations. Associations with false discovery rate (FDR) adjusted at 10% was used to determine the P-values thresholds.

Putative Candidate Gene Analysis and Expression Data

To find candidate genes or putative related proteins of SNP flanking-regions, BLASTx search was conducted for significant MTAs against recently released genome sequence IWGSC RefSeq v1.0⁹. Each MTA was searched for IWGSC sequence information in Ensembl plant for T. aestivum¹⁰. The flanking sequence available for the SNP marker with maximum bases (1,000 bases before and after the SNP) was considered for BLASTx analysis. We also looked at the number of high confidence genes adjacent to the significant MTAs using the RefSeq v1.0.Gene Ontology (GO) annotation of the potential candidate genes was carried out using Blast2GO pro tool v.3.1.3 (Conesa and Götz, 2008). The expression profile of all the putative candidate genes associated with the identified SNPs were checked using wheat RNA-seq expression database of polyploid wheat¹¹. This database consists of the transcript profile of five tissues (grain, leaf, root, spike, and stem) at 3 different time points (growth stages) each and environmental treatments (Pearce et al., 2015). Expression of the gene was measured in units of FPKM (Fragments Per Kilobase of transcript per Million mapped reads). Expression profile was carried out to further provide supporting evidence to corroborate candidate genes (tissue and stage of expression).

Results

Phenotypic Variation and Correlation Analysis

The frequency distribution of phenotypic data of 404 genotypes characterized for 30 traits is given in the Supplementary Figure S1. The phenotypic variations of six quantitative traits (DTH, DTM, PH, SL, Awn_L, and Leaf_L) was recorded in multiple environments. Phenotypic variation of these traits among genotypes was corroborated by mean, standard deviation, range and coefficient of variation (Supplementary Table S3). The mean value of DTH, DTM, PH, SL, Awn_L, Leaf_L varied from 66.72 to 92.0 days, 105.28 to 135.37 days, 94.45 to 113.18 cm, 10.33 to 13.30 cm, 0 to 19.50 cm, 16.25 to 41.40 cm, respectively. This data revealed extensive variation in the traits of the diverse set suggesting the suitability of genotypic panel for association studies. Phenotypic values for each of the six traits were found normally distributed (Supplementary Figure S2). Analysis of variance (ANOVA) was conducted to test the effects of genotype (G), environment (E) and their interactions (G × E). Significant differences were observed among the genotypes (p < 0.0001), the effect of environment and their interaction (G × E) indicating the environmental effect on these traits (Supplementary Table S4). Estimates of correlation coefficients of this combined analysis are shown in Supplementary Table S5 and in Figure 1, a positive correlation was observed for DTH with DTM (0.36), SL (0.18), and PH (0.17) while SL exhibited negative correlation with PH (-0.17).

FIGURE 1

Figure 1. Correlation between six quantitative traits viz. Days to heading (DTH), days to maturity (DTM), spike length (SL), plant height (PH), leaf length (Leaf_L) and Awn length (Awn_L) for 260 genotypes. Significant correlations are designated with asterisk (^∗).

SNP Markers Statistics

Quality preprocessing of 35,143 markers obtained from 35K chip was done by using PLINK software¹² (Purcell et al., 2007). 6,041 monomorphic markers were excluded from the analysis. Out of 29,102 SNP markers, 8,673 SNPs failed frequency test (MAF <0.05) and 1,383 markers removed failing missingness test >0.05. Only 2 individuals for low genotyping (MIND >0.2) were removed. Further, 4,740 SNPs were excluded for lack of their physical position and 146 being duplicate markers. After filtering, 402 genotypes with 14,160 SNP markers were used for GWAS. These markers covered a genetic distance of 4364.79 cM, with an average density of 0.3 cM. Marker density was found highest for B genome (1029.6 markers per chromosome) followed by A (788.9 markers per chromosome) and D genome (207.8 markers per chromosome). Among the genome, chromosome 2B had the highest number of markers (1324) while 4D chromosome spanned the lowest number of markers (55) (Supplementary Table S6).

Population Structure and Linkage Disequilibrium

The mean GD and the PIC for the whole genome were 0.36 and 0.29, respectively. Both GD and PIC of the A genome (0.357 and 0.286) and B genome (0.372 and 0.291) were higher than the D genome (0.345 and 0.276). The number of markers, map length, GD and PIC for each chromosome are shown in Supplementary Table S6.

In the present study, the population structure of a diverse panel of 402 wheat genotypes was investigated on the basis of a ΔK method of model-based Bayesian clustering using 14,160 SNP markers. Population structure analysis clearly indicated the existence of three distinct major subpopulations in the bread wheat panel, which was found consistent with the results of the PCA and neighbor-joining (NJ) tree analysis (Figure 2). Subgroup I, the largest group with 169 accessions, was dominated by recently released varieties and breeding lines adapted to Northern wheat growing zone of the country and genetic stocks for biotic resistance (Rust and Karnal bunt). DPW621-50 (2011), HD2967 (2013), WH1105 (2013), HD3059 (2013), DBW88 (2014), HD3086 (2014), and DBW90 (2014) are some recently released varieties. The pedigree showed that the varieties DPW621-50, DBW88, and HD3059 had common pedigree (Supplementary Table S2). Breeding lines for instance HUW675, HUW666, HPW373, HD3133 and varieties MP1201, HS507, HS542 and WH1105 had MILAN in the parentage. Subgroup II consisted of 87 accessions, mainly comprising local landraces from pre green revolution era; Subgroup III had 146 accessions, predominantly from the warmer region of the country and also comprised of early maturing genotypes (short maturity duration of about 120 days) released for late sowing (toward end of November and to mid of December) in different agro climatic zones viz., K8962, Raj3765, DBW16, and MP3336 having HD2160 (a triple dwarf genotype) as a common progenitor in their background (Supplementary Table S2). Early Mexican cultivars that paved the way for green revolution, Sonalika, SONORA64, Safed Lerma appeared in this cluster along with the derivatives of Sonalika like UP262, HW2001, and Lok54.

FIGURE 2

Figure 2. Population structure and diversity analysis of current GWA study panel. (A) Population structure based on STRUCTURE when K = 3. (B) Neighbor-joining based clustering observed in the study panel using 14160 SNP markers. (C) Three-dimensional plot of the first three principal components, and (D) heat map of pairwise kinship matrix of 402 genotypes.

Linkage disequilibrium (LD) decay distance in the selected panel was found highest in the D genome which decayed at about 5 cM (r² = 0.02) as compared to ∼2 cM in A and B genomes (Supplementary Figure S3). Faster LD decay in D genome vis-à-vis A or B genome has been reported earlier in GWAS of wheat (Lopes et al., 2013; Zhang et al., 2013). With an increase of the genetic distance, the r² value of the A, B, and D genomes decreased gradually. Genome A (62.7%) showed the highest frequency of physically linked locus pairs followed by B (58.0%) and D (53.6%) genomes.

Genome–Wide Association Analysis

In order to detect the most significant marker-trait associations, CMLM was employed to deal with the confounding effect of the population structure. This was followed by the inspection of Q–Q plots and Manhattan plots for evidence of P-value inflation (Supplementary Figures S4, S5). Based on the stringent criterion of -log₁₀ P > 4, we detected 99 significant MTAs ranging from 7.49 e-05 to 2.47 e-11 for 17 qualitative traits (Figure 3, Table 1, Supplementary Figure S4, and Supplementary Table S7) explaining 5.3–33.3% phenotypic variations. It is imperative to note that, not every gene is likely to be represented by 35K SNP array based markers. Therefore the markers in linkage disequilibrium indicates either it is the causative gene itself or might be in close linkage to the causative gene. For color related traits, a total of 22 SNPs were found associated with coleoptile color on five chromosomes, i.e., 2A, 4B, 6B, 5B, and 6A. However, the genomic region on chromosome 6B was represented by eighteen SNPs, mapped within genetic distance of 62.83–67.99 cM (distance interval of 5.15 cM) which collectively explained 23.9% of the phenotypic variation (Supplementary Figure S4 and Supplementary Table S7). For awn and ear color, a significant MTA was detected on chromosome 1B but at different loci, i.e., at 24.64 cM accounting for a phenotypic variation of 8.8% and at 8.24 cM explaining 10.5% phenotypic variation, respectively (Supplementary Figure S4 and Supplementary Table S7).

FIGURE 3

Figure 3. Five different categories of phenol coloration viz. absent, weak, medium, strong, and very strong are shown. These variations in coloration are based on the phenol oxidase activity present in the seed coat. Manhattan plot and Q–Q plot for the trait of grain coloration with phenol (Grn_Ph) as observed in the study.

TABLE 1

Table 1. Genome wide significant MTAs identified for the qualitative traits using CMLM.

For waxiness characters, two MTAs were identified for leaf sheath wax on chromosomes 6D (7.48 e-05), 3A (P < 0.0001) hitherto not reported (Supplementary Figure S4) and one MTA for peduncle wax on chromosome 3A (2.87 e-05). MTAs associated with leaf sheath wax and peduncle wax contributed to the trait negatively. For glume related traits, six MTAs were detected for glume pubescence on chromosomes 1A (2), 1B (1), and 2B (3) explaining phenotypic variation ranging from 8.9 to 12.0% with positive effect. For shoulder width and brush length, significant associations were detected on chromosome 2D and 3A, respectively (Supplementary Figure S4 and Supplementary Table S7). For awn related traits, a genomic region found associated with awn attitude represented by three SNP markers (AX-94613491, AX-94519690, and AX-94453668) on chromosome 5A spanning a region from 59.99 to 70.36 cM. The phenotypic variation explained by SNPs ranged from 24.2 to 25.0% and all the three SNPs showed positive effect on the awn attitude. For awn presence, several markers or regions were identified across the chromosomes (1A, 1B, 1D, 2A, 2B, 3A, 3B, 4A, 4B, 5A, 5B, 6B, 6D, 7A, and 7B) explaining phenotypic variation ranging from 28.7 to 33.9%. A chromosomal region of 12.22 cM (66.99–72.22 cM) on chromosome 5A harbored significant MTAs associated with multiple traits (awn length, auricle color, and awn presence).

For plant growth habit two SNPs were detected on chromosome 1B but at different loci, one at 8.24 cM and other at 38.86 cM indicating the role of two independent loci on chromosome 1B explaining phenotypic variation ranging from 20.0 to 21.2%.

For grain related traits, in the current study 9 SNPs on chromosome 2A spanning 0.71 cM region (124.18–124.89 cM) were found significantly associated with phenol color indicating the importance of this region. These SNPs explained 15.1–18.6% phenotypic variation. For grain texture (phenotype scored as hard or soft), a total of two regions were detected; one on chromosome 7A containing three markers (100.09 cM) and other on chromosome 6B comprising four markers (62.83 cM). Both these MTAs for grain texture contributed negatively to the trait. For germ width, significant association was detected on chromosome 4A. Only one MTA (AX-94670534) for ear shape was detected on chromosome 7B explaining 8.3% of phenotypic variance having negative effect. We did not find any significant MTA for ear density.

For the six quantitative traits (DTH, DTM, PH, Leaf_L, SL, and Awn_L), a total of 47 significant SNPs were identified in five environments which explained 5.2–47.3% of phenotypic variation (Table 2 and Supplementary Figure S5). We successfully detected both previously reported genomic regions and novel loci for the traits in wheat (Tables 1, 2 and Supplementary Figures S4, S5). Flowering time or DTH is a crucial trait which affects the adaptation of wheat in its target environment. A total of 5 SNPs for DTH were detected on chromosome 4A, 5A, and 7D with phenotypic contributions ranging from 19.1 to 32.5%. Out of the 5 SNPs associated with DTH, SNP AX-9454244 on chromosome 4A (78.09 cM) and SNP AX-95187165 (89.02 cM) on chromosome 5A showed pleiotropic effect on DTM. For DTM, a locus on chromosome 2A at 179.61 cM has been detected in the three environments (E2, E3, and E5; Table 2) explaining the average phenotypic variation of 22.9% suggesting the importance of this region while another locus on chromosome 2A at 83.23 cM was observed for three traits (DTM, PH, and SL). Two MTAs were detected on chromosome 5A (AX-95187165, 89.02 cM and AX-95652310, 72.22 cM), indicating the presence of two independent loci on chromosome 5A with a positive effect on DTM.

TABLE 2

Table 2. Significant MTAs identified for quantitative traits in the current and previous study.

For PH, a total of six MTAs were identified, one each on chromosome 1B (24.64 cM), 2A (83.23 cM), 2B (104.59 cM), 7A (29.9 cM), and two on 5D (1.58 cM) considering all the environments (Table 2). MTAs significantly associated with SL were mainly distributed on chromosome 1B, 2A, 3A, 3B, 3D, 5A, 7A, and 7B. The phenotypic variation contributed by SNPs ranging from 8.3 to 27.6% for SL. The SNP AX-94517196 (83.23 cM) on chromosome 2A (Supplementary Figure S5 and Table 2) showed 15.6% of phenotypic variation having positive effect on the trait. Three MTA for flag leaf length were detected on chromosome 4A, 5A, and 7B explaining phenotypic variation ranging from 9.7 to 15.7% (Table 2). The SNP AX-95196340 and AX-94406861 on chromosome 7B and 5A, respectively, showed positive effect on leaf length. For awn length, a total of 8 significant MTAs were identified mainly distributed on chromosome 1B (1), 3B (1), 4B (2), 6D (2), and 7B (2) which collectively explained 18.0% of the phenotypic variation. Two loci, one at 1.72 cM on chromosome 7B and other at 17.65 cM on chromosome 6D shared association with awn length and awn presence.

Pleiotropy Effect

In the present study, we observed same SNPs with multiple traits which could be due to pleiotropy or different causal genes in LD for instance SNP AX-94656878 at 83.23 cM (chromosome 2A) explained variation for two traits (PH and DTM) (Table 2). Similarly, another locus (SNP AX-94527988) on chromosome 3A at 115.9 cM was found pleiotropic with LS_Wax and Ped_Wax (Table 1). Also, the two SNPs associated with DTH on chromosome 4A (78.09 cM) and 5A (89.02 cM) were found linked with DTM (Table 2). The pleiotropic effects observed in the study were in agreement with the Pearson’s correlations observed between the agronomic traits (Figure 1).

Identification of Putative Candidate Genes and Expression Analysis

We identified several putative candidate genes such as storage protein activator (spa), beta-amylase 2 (bmy2), cytochrome P450, shikimate kinase, b-ZIP transcription factor, for the phenotypic variations of the traits (Table 3). These putative proteins identified were highly homologous to different species of Triticum or Aegilops. Highest number of putative candidate genes were observed for MTAs associated with SL encoding a total of five candidate genes [actin-related protein subunit 3 (ARPC3), DIMINUTO, replication protein A, carboxypeptidase D, and basic region/leucine zipper].

TABLE 3

Table 3. SNPs significantly associated with agro-morphological traits and putative candidate genes identified in the study.

To determine the relative expression profile of the identified transcripts in broad range of tissues from different developmental stages, the published RNA-seq data and the Wheat-Exp web tool of the wheat cultivar, Chinese Spring was explored (Choulet et al., 2014; Pearce et al., 2015). The expression profile of significant SNPs encoded putative candidate genes is given in the Table 3 and Supplementary Table S8. FPKM value > 5 was considered for tissue and developmental specific expression check.

Effect of Favorable Alleles on Agronomic Traits

Early maturing, high yielding wheat genotypes are of immense importance toward increasing the cropping intensity as well as ensuring high input use efficiency particularly for inputs like water, which are going to be scarce. Therefore the present study dissecting important agronomic traits such as DTH, PH, and SL enables utilization of available diversity by exploiting associated markers. SNP alleles which led to decrease in DTH, PH and increase in SL were considered as “favorable alleles” and vice-versa was defined as “unfavorable alleles.” Figure 4 depicted higher frequency of favorable alleles which led to decrease in PH and DTH with phenotypic variation of 17.6 and 9.0%, respectively. Similarly, by increasing the number of favorable alleles, SL increased with R² of 6.4%. Results of the study showed that favorable alleles exhibited significant positive effects on the phenotypic traits as compared to the unfavorable alleles. This would help in cultivar adaptation and finally to grain yield.

FIGURE 4

Figure 4. Linear regression between number of (A) favorable alleles (B) unfavorable alleles for spike length (SL), plant height (PH) and days to heading (DTH). Higher regression value (R²) was observed in case of PH when compared to SL and DTH for both favorable and unfavorable alleles.

Traits Sharing Co-localized Genomic Regions

In the present study, the most promising co-localized genomic region was identified on chromosome 1B at 24.64 cM associated with four traits (Glu_Pub, Awn_Col, SL, and PH) and 26.22 cM with Awn_L. The genomic co-location of loci (24.64 cM) with four traits implies either a strong physical linkage between genes underlying these important traits, or a pleiotropic effect. Therefore, to dissect the genetic cause of the observed association, LD patterns and candidate genes underlying the region and transcript profile of the targeted region were investigated (Figure 5). The SL and PH at this locus showed greater LD estimates (>0.8) indicating closely dependent biological processes. Whereas comparatively moderate and low LD was observed with Glu_Pub and Awn_Col, respectively, might be due to low overall recombination vis-à-vis greater recombination frequency with other genomic regions. Notably, this locus harbored three candidate genes for SL. SNP AX-94981940 (-log₁₀ P = 4.21), was annotated as a replication protein A subunit and its transcript expression was almost solely abundant in the young spike at Zadoks 32 stage (Table 3). Replication protein A has important role as single strand DNA binding protein in various DNA metabolic pathways (Aklilu and Culligan, 2016). Similarly, the other SNP encoded a protein carboxypeptidase-D which functions as a positive regulator of grain size in rice (Li et al., 2011). The sequence of SNP AX-94561972 linked with SL annotated as basic region/leucine zipper protein. The bZIP transcription factor family plays an important role in growth, development, and response to abiotic or biotic stresses (Yin et al., 2017). It is interesting to note that PH shared common significant loci with SL showing high correlation between these traits in concurrence with the previous results (Sukumaran et al., 2014). SNP AX-94626335 (-log₁₀ P = 5.10), associated for Glu_Pub at this locus (24.64 cM) was annotated as metal tolerance protein (MTP) which is known for its potential involvement in providing a sink for trace element storage in wheat grains (Vatansever et al., 2017). Earlier Echeverry-Solarte et al. (2015) also reported the influence of glume pubescence on SL by identification of a cluster of co-localizing QTL on same locus for both the traits. Another Glu_Pub associated SNP AX-94664731 on chromosome 1B at 24.64 cM annotated tetratricopeptide repeat protein SKI3, showed highest expression at Zadok39 growth stage of spike. The TaFlo2-A1 gene, an orthologous of rice Flo2 has four motifs of tetratricopeptides found associated with thousand grain weight (Sajjad et al., 2017) and F-box protein containing domains of tetratricopeptides known to regulate plant development and their abundance during spike development in wheat (Hong et al., 2012). The expression patterns of the putative candidate genes in different organs (Figure 5) are consistent with the RNA-seq FPKM expression patterns. The single genomic locus identified for these important related traits, needs further studies to fine map and validate the identity of the causal locus.

FIGURE 5

Figure 5. (A) Genomic region (24.64 cM) associated with Glu_Pub, Awn_Col, SL, and PH on chromosome 1B as identified by GWAS. (B) Pairwise LD estimates in the haplotype block for markers associated with these traits. (C) Graphical representation of a chromosome stretch showing putative candidate genes on a locus 24.64 cM (chromosome 1B). Five candidate genes are represented from 1 to 5 [(1) Basic region/leucine zipper protein, (2) replication protein A subunit, (3) carboxypeptidase D, (4) metal tolerance protein, (5) tetratricopeptide repeat protein SKI3]. (D) Tissue and stage specific expression profile of candidate genes 2, 4, and 5 (Wheat-Exp database).

Discussion

The diversity panel selected in this study has high GD (0.363) and PIC (0.29) indicating higher polymorphism than listed in the previous reports (Liu J. et al., 2017; Eltaher et al., 2018). Further, the B genome had higher GD and PIC followed by A and D genome, consistent with the previous report (Ain et al., 2015). The highest LD decay rate of 5 cM for D genome obtained in this study employing 14,160 SNP markers was found in congruence with 90K SNP (Sukumaran et al., 2014) and 9K (Lopes et al., 2013) marker data. The results from the three clustering methods (Structure, PCA, and NJ tree analysis) showed the presence of three subpopulations in this study consistent with the geographic origins and pedigrees of the selected panel. Thus the marker density, diversity and sample size of this study is sufficiently powered to capture allelic variations for the selected traits. Ma et al. (2013) reported that various imputation methods could be used to impute the data from low density to high density, i.e., from 3K to 54K, and subsequently from 54K to 777K. Therefore the data generated with 35K breeders array can be imputed to high density using 820K information.

In this study, a GWAS panel was characterized for 36 agro-morphological traits identified 146 MTAs (-log₁₀ P ≥ 4) for 23 traits. For majority of the heritable traits, at high significant level single locus has been identified indicating that they are controlled by small number of loci, for instance PGH (1B), LS_wax (6D, 3A), Ped_Wax (3A), Sh_Wid (2D), Ear_Col (1B), Ear_Sh (7B), Awn_Col (1B), Awn_Att (5A), Brush_L (3A), Germ_Wid (4A), and Grn_Ph (2A). For coleoptile color, the genes that regulate anthocyanin biosynthesis pathway, have been cloned and mapped on homoeologous groups 3 and 6 (phenylalanine ammonia-lyase), homoeologous groups 1 and 2 (chalcone synthase), homoeologous group 5 (chalcone-flavanone isomerase) (Li et al., 1999). The presence of the dominant allele at the Rc-1 homeologous loci responsible for anthocyanin pigmentation in coleoptile was correlated with F3H (flavanone 3-hydroxylase) gene on chromosome 2A (Khlestkina et al., 2008). However, in this study, besides chromosome 2A, we also detected the loci for coleoptile color on chromosome 6B, 5B, 4B, and 6A. In agreement to this, Sutka (1977) identified the gene designated Rc4 for coleoptile color on chromosome 6B, however, it was not further confirmed in any study (Khlestkina et al., 2002; McIntosh et al., 2014) while suppressors playing role in the intensity of the coleoptile coloration were identified on chromosomes 2A, 2B, 2D, 4B, and 6A of “Mironovskaya 808” bread wheat variety. Hence the loci identified herein further confirmed the role of chromosome 6B for coleoptile color. For Awn_Col and Ear_Col, we identified loci on chromosome 1B but at different positions, i.e., 24.64 and 8.24 cM, respectively. Earlier, Zeven (1983) also reported a semi-dominant gene (Rg) on chromosome 1B responsible for the brown ear character of bread wheat. For Aur_Col, contrasting to the region reported on chromosome 4A and 5B (Zeven, 1985), we detected its locus on chromosomes 4B and 5A owing to the instability of its expression.

Epicuticular wax is associated with increased drought tolerance in wheat (Bennett and Schnurbusch, 2016), rice (Haque et al., 1992), maize (Meeks et al., 2012), barley (Febrero et al., 1998), and many other crops. Herein, we report an additional locus for LS_Wax on chromosome 6D besides previously identified genomic region on chromosome 3A for waxiness. Interestingly, the SNP AX-94527988 (chromosome 3A) was found associated with both LS_wax and Ped_wax thereby indicating its pleiotropic behavior. In fact, this result indicated that some casual gene(s) might exist in this genomic region for wax, as the common MTA AX-94527988 annotated cytochrome P450 protein which leads to a double-hydroxylation to the corresponding oxo-2-alkanol esters which are also previously detected for both peduncle and flag leaf waxes (Racovita et al., 2016).

Glume pubescence appears to have a beneficial influence on drought/cold tolerance (Borner et al., 2005). In this study, we identified significant MTA (1.47 e-07) on chromosome 1A which is in agreement with Sears (1953) who mapped a gene (Hg) responsible for Glu_Pub on chromosomes 1A. In addition to this, we detected genomic regions associated with Glu_Pub on chromosomes 1B (24.64 cM) and 2B (76.24, 76.38, and 104.59 cM) which might be considered as novel region controlling the trait. MTA AX-95023665 linked to Glu_Pub encoded fatty acid biosynthetic process. Glas et al. (2012) observed that methyl ketones which are produced during fatty acid biosynthesis were the major constituent of type VI trichomes of the wild tomato Solanum habrochaites f. glabratum and are very effective in protecting the plant from pests. Pubescent plants also produced a higher number of grains per spikelet compared to non-pubescent plants (Maes et al., 2001).

Several MTAs detected for awn presence were found distributed across the wheat chromosomes except for 2D, 3D, 4D, 5D, 6A, and 7D. However, we could not locate any awn development dominant inhibitor genes Hd, B1, B2 fine mapped on chromosome 4AS, 5A, and 6BL in hexaploid wheat (Sourdille et al., 2002; Yoshioka et al., 2017) which may be due to skewness for the awned genotypes in the diversity panel. For awn attitude, one QTL located on chromosome 5A (59.99–70.36 cM) was identified. To the best of our knowledge, none of the previous studies have reported a genomic region for Awn_Att, suggesting this could be a responsible locus for the trait. The MTA AX-94613491 encoded hexose carrier protein HEX 6 which is responsible for controlling the flux of carbon and plays a role in the carbohydrate transport and distribution in plant cells¹³.

A gene for high PPO activity responsible for grain color was mapped on the long arm of chromosome 2A in wheat mapping population (Simeone et al., 2002). Similarly, in the current study, 9 SNPs on chromosome 2A spanning 0.71 cM region (124.18–124.89 cM) were found significantly associated with phenol color indicating the importance of this region. These SNPs explained 15.1–18.6% phenotypic variation. Out of 9 SNPs associated with phenol color, 6 SNPs encoded 5 proteins (Table 3). The transcript profile of nucleoredoxin 1-2 (NRX1) gene encoded by SNP AX-94738314 associated with phenol color showed highest expression (FPKM-272.53) in grain at Zadok 72 stage (early milk). Protective effect of NRX1 boosted the H₂O₂ detoxification capacity of catalase, thereby protecting the plant cell from oxidative stress (Kneeshaw et al., 2017). Significant SNPs associated with grain texture (phenotype scored as hard or soft), encoded four proteins (Table 3). Most of the proteins identified in this study for grain traits were similar to proteins reported by Arora et al. (2017) for grain related traits.

The genetic architecture of the quantitative traits is complex as controlled by many loci with small effect. Several significant markers mined for six complex traits in this study were co-localized with the previously reported QTL regions (Table 2). Flowering is controlled by a complex network of genes integrating vernalization response genes (Vrn) on chromosomes 5A (Vrn-A1and Vrn-A2) and 7BS (Vrn-A3), photoperiod response gene on chromosome 2, and earliness genes on chromosomes 1A and 3A (Fowler and Laudencia-chingcuanco, 2016). The MTAs for DTH were detected on three chromosomes 4A, 5A, and 7D. MTA AX-95187165 (89.02 cM) identified for DTH on chromosome 5A was mapped at nearly same position of the functional gene Vrn-A1 (90 cM) indicating its role for determining DTH (Sukumaran et al., 2014). Another putative candidate, AX-94542441 (chromosome 4A) associated with DTH is an ortholog of shikimate kinase-like protein. Likewise in rice (Oryza sativa) two shikimate kinase isoforms OsSK1 and OsSK3 accumulate to high levels during the heading stage of panicle development and involved in floral organ development (Kasai et al., 2005). The most significant association for DTH was found on chromosome 7D stable in four environments (E1–E4). On similar lines, Lopes et al. (2013) also reported a significant QTL for DTH on chromosome 7D, associated for more than 30% of DTH variation but at different position.

The SNPs significantly associated with DTM were identified on chromosomes 3B, 5A, 5B, and 6A corresponding to the earlier reported genomic regions for DTM on chromosome 3B (Sukumaran et al., 2014, 2018; Okechukwu, 2017), 5A (Gahlaut et al., 2017), 5B (Gahlaut et al., 2017; Zou et al., 2017), and 6A (Sukumaran et al., 2014). The significant QTL for DTM harboring three SNPs were observed on chromosome 3B (61.38, 84.9, and 85.27 cM) were closely co-localized with the previously identified MTAs for DTM, indicating that these QTLs were stable and could be detected in different environments. Another noteworthy region, associated with DTM was identified on chromosome 2A (179.61 cM) consistent in three environments encoding BTB/POZ domain and ankyrin repeat-containing protein which plays key role in plant growth and development stages (Sharma and Pandey, 2016). Moreover, chromosome 2A associated with DTM also encompasses a region (83.23 cM) governing multiple traits (PH and SL) thereby representing the correlation of these traits in diverse panel in agreement with the previous reports (Kamani et al., 2017; Zhai et al., 2016).

For PH, as many as six MTAs were identified on chromosome 1B, 2A, 2B, 5D (2), and 7A. The genomic region of MTA (AX-94941145) identified at 29.9 cM on chromosome 7A falls in the region of the reduced height gene Rht22 (Xgwm471-29.5 cM, Xgwm350-20.1 cM) (Peng et al., 2011). Similarly, the MTA identified for PH on chromosome 2B (104.59 cM) and 5D (1.58 cM) found in proximity to the reduced height genes Rht4 (Xwmc317-106 cM) (Ellis et al., 2005) and Rht23 (Xgdm63-4.7 cM, Xbarc110-11.1 cM) (Chen et al., 2015). We also observed an MTA on chromosome 2A where Rht7 gene had already been reported but at different map position. For two correlated traits PH and DTM (Mohibullah et al., 2012), a significant MTA (AX-94656878) was detected on chromosome 2A annotated bZIP transcription factor. In plants, these factors regulate genes in response to abiotic stress, seed maturation, flower development and pathogen defense (Jakoby et al., 2002). Similarly, several studies reported a moderate, but significant correlation between heading time and PH (Sultana et al., 2002; Mohibullah et al., 2012). The MTA AX-94941145 identified on chromosome 7A annotated probable LRR receptor-like serine/threonine-protein kinase At3g47570. Further, we investigated the possibilities of semi dwarfing genes on chromosome 4B and 4D which is present in Indian cultivars (Sheoran et al., 2013) but could not detect any MTA linked to these genes (Rht-B1b and Rht-D1b) suggesting either these genes were eliminated during filtering or may not reach significant threshold level. Similarly, Ain et al. (2015) did not find any MTA for semi dwarfing genes on chromosome 4B and 4D employing iSelect 90K SNP chip. Several MTAs for PH have been reported previously on chromosomes 2A (Ain et al., 2015; Mengistu et al., 2016), 2B (Zanke et al., 2014; Ain et al., 2015; Gao et al., 2015) and 7A (Gao et al., 2015;Soriano et al., 2017).

Wheat domestication genes Q, compactum (C), sphaerococcum (S1) related to spike morphology have been detected on chromosomes 5A, 2D, and 3D, respectively (Johnson et al., 2008). In the present study, eight loci were detected on chromosomes 1B, 2A, 3A, 3B, 3D, 5A, 7A, and 7B for SL which were partially consistent with those of Zhou et al. (2017), who reported QTLs for SL on chromosome 3A, 3B, 5A and with those of Ma et al. (2007), who reported genomic regions for SL on chromosome 1B and 7A. These results indicated that multiple loci having unequal effects can influence the variations in the SL. It is interesting to note that PH shared common significant loci (83.23 cM) with SL showing a high correlation between these traits in concurrence with the previous results (Sukumaran et al., 2014). Besides MTAs on chromosome 1B (discussed earlier), another SL associated SNP AX-94722223 (chromosome 5A) harboring actin related protein (ARPC3) which has been known to play a key regulator of cytoskeleton dynamics controlling multiple developmental processes in a variety of tissues and cell types (Qi et al., 2017). We expected genes contributing to variation in SL to be most strongly expressed within the different growth stages of spike. In fact, 5 putative candidate genes identified for SL showed expression FPKM >5 in spike tissue depicting highest expression (FPKM 143.65) at growth stage 65 (Zadok’s Scale), corroborating its causal effect.

Awns were reported superior to the flag leaves on a cellular and physiological level during the grain filling period contributing 40–80% of the photosynthetic assimilates accumulated in the wheat grain (Li et al., 2006). Most significant MTAs for Awn_L were reported on several chromosomes viz., 1A, 1B, 2B, 2D, 3B, 4A, 5B, 6D, and 7B (Wu et al., 2016; Liu Y. et al., 2018). However, some of the chromosomal regions associated with Awn_L for instance 1.72 cM (7B) and 85.27 cM (3B) were detected for both Awn_L and SL. This study further corroborated the result of Echeverry-Solarte et al. (2015), who reported a novel QTL for Awn_L on chromosome 7B harboring two consistent loci associated to supernumerary spikelet (SS) and putative QTLs for PH, DTH, DTM making it an important loci for future studies. In the present investigation, cell elongation protein DIMINUTO predicted for Awn_L associated SNP locus (AX-95025537 and AX-95012310) on chromosome 7B (Supplementary Table S8) have been implicated in regulating cell elongation (Takahashi et al., 1995). The genomic regions that contributed to Leaf_L were found associated with chromosomes 4A and 5A with phenotypic variations ranging from 15.5 to 15.7%. In earlier reports several chromosomal regions viz., 2B, 3A, 4A, 4B, and 5A were detected for flag Leaf_L (Jia et al., 2013; Wu et al., 2016). Positive and significant correlation between flag Leaf_L and SL revealed their role in increasing yield (Wu et al., 2016).

Mining of superior/favorable alleles is essential for improving the complicated earliness trait in wheat using marker assisted selection. In recent years, association mapping has been widely used in exploring the elite alleles of many agronomic traits such as yield related traits (Sun et al., 2017), heading days and PH (Ain et al., 2015; Ogbonnaya et al., 2017) and water soluble carbohydrates (Dong et al., 2016) in wheat. In the present study, the phenotypic effect value of the favorable alleles of DTH, PH, and SL was evaluated and inferred to have positive effect on the respective traits. The candidate genes and the SNPs linked with the economic important traits identified in this study could help in designing new strategies to hoard superior alleles for these key traits in future marker based breeding. Some novel regions identified in the present investigation could have been previously detected, but comparison of the positions of the SNPs linked to the respective traits was not possible due to the limitations of the various marker system used in different studies.

This study identified 146 MTAs for 23 agro-morphological traits (Supplementary Table S9), and putative candidate genes using the recently released genome sequence by IWGSC RefSeq v1.0 (Appels et al., 2018). MTAs specific to less explored traits such as awn length and glume pubescence were targeted for visualizing future needs of breeders in developing efficient and resilient wheat varieties. The chromosomal region controlling multiple traits were also identified which should pave the way for selection and may prove effective for pyramiding favorable alleles. Here we discovered novel candidate genomic regions together with previously reported genes which require further validation and testing in the wheat germplasm. Therefore, the significant MTAs identified having known candidate genes are being subjected to conversion as Kompetitive Allele Specific PCR (KASP) markers that can be efficiently used to transfer alleles into elite wheat genotypes (Rasheed et al., 2016). These useful genomic resources and PCR based markers (KASP markers) could be utilized for introgression of traits through marker assisted selection (MAS). These will strongly enhance systematic study of the genetics, comparative genomics and evolution of wheat, and will expedite isolation and characterization of genes controlling agronomically important traits, such as yield and abiotic stress.

Author Contributions

RT and DiK conceived the theme of the study. SS, DeK, NR, RuS, SusP, SJ, MI, AJ, NK, UA, and SurP did the computational analysis. SS, SJ, MI, RaS, PS, RT, and DiK drafted the manuscript. DeK, NR, RuS, CS, and AG did the phenotyping. NR and DeK contributed in wet lab work. SJ, MI, DiK, AR, GS, and RT edited the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by Indian Council of Agricultural Research; Ministry of Agriculture and Farmer’s welfare, Government of India by providing financial assistance in the form of CABin grant (F. No. Agril. Edn.4-1/2013-A&P) to ICAR-IASRI and ICAR-IIWBR.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors are thankful to Indian Council of Agricultural Research, Ministry of Agriculture and Farmer’s welfare, Government of India for Advanced Supercomputing Hub for Omics Knowledge in Agriculture (ASHOKA) facility at ICAR-IASRI, New Delhi, India created under National Agriculture Innovation Project, funded by World Bank.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.00527/full#supplementary-material

Footnotes

References

Ain, Q., Rasheed, A., Anwar, A., Mahmood, T., Imtiaz, M., Mahmood, T., et al. (2015). Genome-wide association for grain yield under rainfed conditions in historical wheat cultivars from Pakistan. Front. Plant Sci. 6:743. doi: 10.3389/fpls.2015.00743

PubMed Abstract | CrossRef Full Text | Google Scholar

Aklilu, B. B., and Culligan, K. M. (2016). Molecular evolution and functional diversification of replication protein A1 in plants. Front. Plant Sci. 7:33. doi: 10.3389/fpls.2016.00033

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, A. M., Winfield, M. O., Burridge, A. J., Downie, R. C., Benbow, H. R., Barker, G. L. A., et al. (2017). Characterization of a wheat breeders’ array suitable for high-throughput SNP genotyping of global accessions of hexaploid bread wheat (Triticum aestivum). Plant Biotechnol. J. 15, 390–401. doi: 10.1111/pbi.12635

PubMed Abstract | CrossRef Full Text | Google Scholar

Appels, R., Eversole, K., Feuillet, C., Keller, B., Rogers, J., Stein, N., et al. (2018). Shifting the limits in wheat research and breeding using a fully annotated reference genome. Science (80-) 361:eaar7191. doi: 10.1126/science.aar7191

PubMed Abstract | CrossRef Full Text | Google Scholar

Arends, D., Prins, P., Jansen, R. C., and Broman, K. W. (2010). R/qtl: high throughput multiple QTL mapping. Bioinformatics 26, 2990–2992. doi: 10.1093/bioinformatics/btq565

PubMed Abstract | CrossRef Full Text | Google Scholar

Arora, S., Singh, N., Kaur, S., Bains, N. S., Uauy, C., Poland, J., et al. (2017). Genome-wide association study of grain architecture in wild wheat Aegilops tauschii. Front. Plant Sci. 8:886. doi: 10.3389/fpls.2017.00886

PubMed Abstract | CrossRef Full Text | Google Scholar

Arruda, M. P., Brown, P. J., Lipka, A. E., Krill, A. M., Thurber, C., and Kolb, F. L. (2015). Genomic selection for predicting head blight resistance in a wheat breeding program. Plant Genome 8, 1–12. doi: 10.3835/plantgenome2015.01.0003