Genetic Dissection of Seedling Root System Architectural Traits in a Diverse Panel of Hexaploid Wheat through Multi-Locus Genome-Wide Association Mapping for Improving Drought Tolerance

Cultivars with efficient root systems play a major role in enhancing resource use efficiency, particularly water absorption, and thus in drought tolerance. In this study, a diverse wheat association panel of 136 wheat accessions including mini core subset was genotyped using Axiom 35k Breeders’ Array to identify genomic regions associated with seedling stage root architecture and shoot traits using multi-locus genome-wide association studies (ML-GWAS). The association panel revealed a wide variation of 1.5- to 50-fold and were grouped into six clusters based on 15 traits. Six different ML-GWAS models revealed 456 significant quantitative trait nucleotides (QTNs) for various traits with phenotypic variance in the range of 0.12–38.60%. Of these, 87 QTNs were repeatedly detected by two or more models and were considered reliable genomic regions for the respective traits. Among these QTNs, eleven were associated with average diameter and nine each for second order lateral root number (SOLRN), root volume (RV) and root length density (RLD). A total of eleven genomic regions were pleiotropic and each controlled two or three traits. Some important candidate genes such as Formin homology 1, Ubiquitin-like domain superfamily and ATP-dependent 6-phosphofructokinase were identified from the associated genomic regions. The genomic regions/genes identified in this study could potentially be targeted for improving root traits and drought tolerance in wheat.


Introduction
Globally, wheat (Triticum aestivum L.) is one of the most important cultivated food grain crops with an area and production of 218.54 million hectares and 771.71 million tons, respectively [1]. Wheat ranks second after rice (Oryza sativa L.) as a major staple food crop in India with an acreage of 30.55 million hectares, and total production and productivity of 107.18 million tons and 3508 kg/ha, respectively [2]. Wheat is adapted to various climatic conditions including areas vulnerable to heat, salt and drought stress [3,4]. Drought is one of the major challenges that limits crop growth and yield in different areas of the world [5,6]. Drought stress can cause significant yield reduction in wheat and its impact varies with intensity, timing and duration of stress relative to crop growth stages. Yield reduction to the tune of 21.0%, 25.8% and 32.0% were reported in wheat under mild (55-70% relative soil water content), moderate (35-55%, relative soil water content) and severe (<35% relative soil water content) drought stress, respectively [7].
Root traits are important not only for water absorption from a drought stress perspective, but also for soil nutrient uptake and environmental stress tolerance [8][9][10]. Morphology of the root system includes basic features such as root length, root length density, root diameter, surface area and volume which influence the root structure's spatial arrangement and are significantly correlated with water uptake and nutrient absorption [11][12][13]. Drought tolerance is directly associated with root diameter, as thicker roots with large xylem vessels are more effective towards extraction of water and nutrients from deep soil layers under rainfed condition [12]. Fine roots and the number of root tips, the major components of root systems, often increase the water and nutrient uptake efficiency by increasing the root surface area and volume [8,10]. Further, deep rooting, an essential root feature allowing access to water from deeper soil profiles, improves crop productivity [14]. Therefore, one approach to improve drought tolerance is to identify and introgress deeper rooting alleles in shallow rooted and drought prone cultivars [15][16][17]. Thus, the study of root system architecture (RSA) features can help define proxy traits to enhance tolerance to various soil types and status of moisture and nutrient stress conditions [18]. RSA characteristics are governed by polygenes with additive cumulative gene action and are influenced by the environment [12]. Due to the difficulty of high-throughput phenotyping of RSA characteristics in field conditions, their optimization was ignored [10,12,19]. However, seedling level root architecture is associated with adult plant root architecture and is being studied to translate the knowledge for later growth stages [20,21].
The knowledge on extent of genetic variability in wheat germplasm collections could pave the way for successful conservation and utilization in crop improvement. In recent years, cost effective advanced high throughput genotyping approaches, such as single nucleotide polymorphism (SNP) arrays and sequencing based genotyping methods such as genotyping by sequencing (GBS), restricted site associated DNA sequencing (RAD-seq) and whole genome re-sequencing, have greatly accelerated genomic level characterization of genetic diversity of crop germplasm resources in many species. Availability of these techniques has made it feasible to conduct whole genome scan studies, such as association mapping (AM) and genomic selection in a variety of crops including wheat, maize (Zea mays L.), and rice [22,23]. Association mapping is an alternative approach to mapping quantitative trait loci (QTLs) controlling complex traits which is based on the principle of linkage disequilibrium between marker and trait. In contrast to the conventional QTL approach, the AM approach enables screening of a large number of alleles present in natural populations of species, thus providing greater resolution mapping of a target trait. In the past couple of years, AM has been widely applied to map traits in many crop species including wheat, rice, barley (Hordeum vulgare L.), maize and soybean (Glycine max L. Merril) [21,[24][25][26][27] with the availability of high-throughput genotyping and phenotyping platforms [22,23].
Exploring genetic variability of root traits in a diverse germplasm collection could assist in development of varieties with desired root features for drought tolerance or target environments. To date, only a few reports are available on genetic variability and QTL mapping of root traits in wheat [28][29][30]. The wheat core and mini core collections constituted from the large wheat accessions conserved in the national and international gene banks could be a valuable genetic resource for mapping complex traits using a genome wide association study (GWAS). The GWAS could facilitate uncovering genomic regions/genes controlling various root architecture traits. In our earlier study, a core set of wheat genotypes was identified from the entire wheat collection in the national gene bank of India based on 35 agro-morphological traits [31]. This core set was subsequently reduced to a mini core collection that captured huge variability for various traits and could be used for mapping root associated traits. The AM analysis on this mini core collection could facilitate identification of genes/genomic regions associated with RSA traits that could be exploited in a breeding program to develop drought tolerant cultivars. Keeping this in mind, the objectives of this research were: (1) to study genetic variation of seedling root system architecture in a diverse set of wheat genotypes including a subset of wheat mini core collection; and (2) to identify genomic regions/candidate regions linked with these traits using association mapping.

Phenotypic Analyses
There were large differences for various root traits (  (Table 1). We observed a high range of coefficient variation (%CV) for all the phenotypic traits. The highest CV was observed for SOLRN (64.10%) followed by FOLRN (44.50%) and LRS (38.39%). The lowest CV was observed for AD (8.46%) followed by SRN (10.29%). The estimates of broad-sense heritability were low to moderately high, ranging from 15.16% for LRD, to 78.95% for AD (Supplementary Table S2).
The frequency distribution pattern of traits is presented in Supplementary Figure S1. Most of the traits followed normal distribution, except SRN and SL which were negatively skewed, while SOLRN and SRL were positively skewed. Selected wheat genotypes with high and low mean performance of root and shoot morphological traits are shown in Supplementary Table S3. The wheat genotype IC28755 exhibited maximum root length of 37.50 cm, whereas the genotype EC6903, a collection from the USA, showed a maximum shoot length of 40.76 cm. The genotype IC128151 (C306) had highest value of 60.41 mg for SDW and the genotype EC576889 had maximum value for RDW at 23.06 mg. The trait root-shoot dry weight ratio indicated the amount of root dry mass to shoot dry mass. The genotype IC532019 had the highest ratio of 0.59 contributing higher root dry mass than shoot at seedling stage. Donors for multiple traits were also identified. Genotype IC82425A exhibited the highest mean values for SOLRN, LRS and RV. Further the genotype IC128151 (C306) was among the top ten genotypes for seven traits, namely TRS, FOLRN, LRD, RLD, RV, SDW and RDW. The genotype IC542076 ranked in the top ten genotypes for seven traits: TRS, SOLRN, RLD, RV, SRN, SDW and RDW, and was therefore considered promising for tolerating drought by extracting water stored in the deeper soil layers. The genotypes IC128151, IC82425A and IC542076 were superior for several root and shoot traits in the drought experiment carried out at the seedling stage under pot screening during the winter season 2019-2020(data not presented). In addition, genotypes with low mean performance were also identified so they could be utilized to generate bi-parental mapping populations or multiparent advanced generation intercross (MAGIC) populations. An exotic collection from the USA, EC339632, was observed to have lowest mean performance for TRS, FOLRN, SOLRN, RLD, RV, SRN, RDW and RL, whereas EC187159 showed lower mean values for TRS, SOLRN, LRS, RLD, SRN, SDW, RDW and SL.   . Correlations between root system architecture traits and shoot traits. TRS, total root size (cm); FOLRN, first order lateral root number; SOLRN, second order lateral root number; LRS, lateral root size; LRD, lateral root density (cm −1 ); RLD, root length density cm −2 ; AD, average diameter (mm); RV, root volume (cm 3 ); SRN, seminal root number; SDW, shoot dry weight (mg); RDW, root dry weight (mg); SL, shoot length (cm); RL, root length (cm); RSDWR, root shoot dry weight ratio; SRL, specific root length (cm mg −1 ). Red boundary for each pairwise correlation displays density ellipses covering 95% points between each variable. The density ellipses are a graphical indicator of the correlation between two variables. It collapses diagonally as the correlation between the two variables approaches either 1 or −1. It is more circular if the two variables are less correlated.

Principal Component Analysis
Principal component analysis (PCA) was performed for 15 RSA and shoot traits to determine their contribution to overall variability. The first four principal components (PCs) in the PCA analysis with eigenvalues >1 contributed a cumulative variance of 79.24% among the genotypes (Supplementary Table S4). The scree plot showing the eigenvalues on the y-axis and the number of components on the x-axis also depicted the number of relevant PCs to explain maximum variance (Supplementary Figure S2A). Based on the contribution towards the principal components, RV, RDW, SOLRN, RLD, TRS, RL, LRS, FOLRN, SRN and RSDWR were the major contributors toward PC1 and were most significantly associated with the genotypic variations. Likewise, SDW and SL were major contributors to PC2, LRD and FOLRN to PC3 and SRL to PC4 (Supplementary Table S4). Length of coordinates or distance from the origin revealed the amount of variation contributed by each trait (Supplementary Figure S2B) which also ascertained about their major contributions. The PCA also indicated trait association, as the correlated traits have fewer angles between them or are mostly parallel to each other. In this case, the positive quadrant had correlated traits SDW and SRN. The results indicated that AD and SRL were opposite to other traits in the second quadrant, hence their negative association with them.

Clustering of Genotypes
The diverse association panel of the wheat mini core subset (136 accessions, four had a poor SNP call and were therefore excluded from the analysis) was divided into six clusters (Supplementary Figure S3) when subjected to hierarchical analysis based on the 15 root and shoot traits, which showed that the wheat genotypes were sufficiently diverse. Among different classes, the largest cluster (II) contained 38 member genotypes and the smallest cluster (III) had 2 members (EC339632, a USDA collection and EC578153 and a synthetic line). Clusters six and three consisted of genotypes with high and low performance respectively for many traits such as TRS, FOLRN and SOLRN (Supplementary Table S5). A heatmap generated by two-way clustering also depicted that cluster six had high performing genotypes, whereas cluster 1 had low performing genotypes for many of the traits (Supplementary Figure S3).

Marker Distribution, Linkage Distribution and Population Structure
The SNP genotyping of 136 wheat genotypes yielded a total of 16,616 polymorphic SNPs with exact chromosomal position and minor allele frequency (MAF) ≥ 10. The SNP markers were proportionally distributed across the 21 chromosomes and their number ranged from 287 for Chr4D to 1179 for Chr2B, and the average SNP markers per chromosome was 791. The polymorphic SNPs covered 14066.28 Mb of wheat genome. The sub-genome A included 5355 SNPs, the sub-genome B had 6343 SNPs and 4918 SNPs were located on the sub-genome D. On sub-genome A, maximum SNPs were found on 2A (952) followed by 1A (911). Further, on sub-genome B, maximum SNPs were found on Chr2B (1178), followed by Chr5B (1125). Chr2D of Sub-genome D had maximum SNPs (1037) followed by Chr1D (895) ( Table 2). The SNP marker density per chromosome per Mb ranged from 0.56 for 4D to 1.59 for 2D, with a mean of 1.182. These polymorphic 16,616 SNP markers were used to estimate genome-wide linkage disequilibrium (LD) and the pairwise pattern of LD was calculated using r 2 of allele pairs between two loci. Among these, a total of 19,528 SNP pairs were in perfect LD (r 2 = 1), with B sub-genome having a maximum (7287) number followed by D (7100) and A (5141) sub-genomes. Overall, in the current association panel, low mean LD was observed across the 21 chromosomes with values ranging from 0.087 (4D) to 0.298 (1D). The average LD (r 2 ) per chromosome ranged from 0.087 for 4D to 0.298 for 1D. The measures of LD decay for all 21 chromosomes are summarized in Supplementary Figure S4. The LD decay in the AM panel was determined at r 2 = 0.157 (background LD). LD decayed the fastest in sub-genome A followed by sub-genome D and B. In sub-genome A, LD decayed at 3 Mb compared to 3.4 Mb in sub-genome B and 4.2 Mb in sub-genome D ( Table 2 and Supplementary Figure S4). Population analysis of the association panel was performed using both STRUCTURE and PCA which indicated optimal number of population (K) as two ( Figure 2). Of the total accessions, 40 accessions belonged to subpopulation 1 and 96 accessions belonged to subpopulation 2. Further, in subpopulation 1 (SP1), 30.50% of the accessions were admixtures and the remaining 62.50% of the genotypes were pure (>80% similarity) which included 25% EC (exotic) collection and 75% IC (indigenous) collection. The IC belonged to multiple states in India (Rajasthan, Gujarat, Haryana, West Bengal, Uttarakhand, Uttar Pradesh, Himachal Pradesh, Madhya Pradesh and Maharashtra). The EC belonged to Israel, China, the USA and Australia. In subpopulation 2 (SP2), 42.70% of the genotypes were admixtures and 57.29% of the genotypes were pure which included 31.25% EC and 68.75% IC collection. The EC mainly belonged to Mexico, Russia, Finland and the USA. Based on PCA analysis, genotypes contributing more variation to the PC1 were parallel to the PC1 axis and genotypes contributing more variation to the PC2 were parallel to the PC2 axis ( Figure 2). A total of 88.2% cumulative variance were explained by two PCs. It also indicated that genotypes used in this study were divided into two structured subpopulations.

Colocalization of Root Architecture Loci
Different root traits are mostly correlated and complex biological mechanisms might be involved in coordination for their expression. The pleiotropic action of genetic loci on different traits or their tight linkage results in a correlation between traits. A total of 11 QTNs were found to be pleiotropic markers (Table 4). A locus on 1D (AX-94448890) was associated with root morphology (SOLRN and RL) and another locus on 2A (AX-94952472) was associated with TRS and RLD. Similarly, other loci were also identified as pleiotropic on 2B (AX-95227366) for RV and RL, 3B (AX-94516395) for SOLRN and RLD, 4A (AX-95105488) for TRS, RLD and RV, 6A (AX-95244609) for SOLRN, LRS, AD and 6B (X-94502864) for RLD, RV and RL. We also identified two pleiotropic loci on the short arm of 7B (AX-95123855) associated with root morphology traits (FOLRN, SOLRN, LRS) and QTN-7B-AX-95119337 associated with TRS and RLD whereas its long arm had a locus (AX-94528392) associated with RL and LRS. Other pleiotropic QTNs were 1DS-AX-94448890 (SOLRN and RL) and 7AL-AX-94664277 (RV and SDW) (Figure 4 and Table 4). These identified loci influencing several traits could be potential markers for marker-assisted breeding after validation.

Colocalization of Root Architecture Loci
Different root traits are mostly correlated and complex biological mechanisms might be involved in coordination for their expression. The pleiotropic action of genetic loci on different traits or their tight linkage results in a correlation between traits. A total of 11 QTNs were found to be pleiotropic markers (Table 4). A locus on 1D (AX-94448890) was associated with root morphology (SOLRN and RL) and another locus on 2A (AX-94952472) was associated with TRS and RLD. Similarly, other loci were also identified as pleiotropic on 2B (AX-95227366) for RV and RL, 3B (AX-94516395) for SOLRN and RLD, 4A (AX-95105488) for TRS, RLD and RV, 6A (AX-95244609) for SOLRN, LRS, AD and 6B (X-94502864) for RLD, RV and RL. We also identified two pleiotropic loci on the short arm of 7B (AX-95123855) associated with root morphology traits (FOLRN, SOLRN, LRS) and QTN-7B-AX-95119337 associated with TRS and RLD whereas its long arm had a locus (AX-94528392) associated with RL and LRS. Other pleiotropic QTNs were 1DS-AX-94448890 (SOLRN and RL) and 7AL-AX-94664277 (RV and SDW) ( Figure 4 and Table 4). These identified loci influencing several traits could be potential markers for marker-assisted breeding after validation.

Allelic Effects of Identified Genomic Regions on Respective Phenotypes
Seven major loci detected in four or more ML-GWAS methods were analyzed to determine the range of phenotypic variations for all traits ( Figure 5). Association panel genotypes were divided into two classes according to allele types. It was observed that all seven QTNs demonstrated a significant effect on respective traits (p ≤ 0.01). Among these seven QTNs, three showed a  . Boxplot for seven reliable QTNs with significant effects (p < 0.01) on corresponding root and shoot traits. For each QTNs, the germplasm lines were divided into two groups according to superior and inferior allele type. The x-axis represents the two alleles for each QTNs, while the y-axis corresponds to phenotypic value. Subfigures, (A-C) represent allelic differences for SOLRN and superior alleles for QTNs Q.SOLRN-1DS, Q.SOLRN-1BL and Q.SOLRN-6AL were T, A and T respectively. Subfigures (D-G) represent allelic differences for LRS, RV, SDW and AD for QTNs Q.LRS-7BS, Q.RV-5BL, Q.SDW-4AL and Q.AD-6AL with superor alleles as T, C, T and C respectively.

Genes Linked to Quantitative Trait Nucleotides
Out of the total SNPs associated with various root and shoot traits, 52 SNP markers were found to be located within genes (Supplementary Table S6). Some important proteins encoded by these genes included glycoside hydrolase superfamily protein, carbonic anhydrase superfamily protein, auxin response factor, nitrogen regulatory protein, transmembrane protein, enhancer of polycomb-like protein and glutamate receptor. More-over, we found that some of the loci corresponding to large effect QTNs and detected by multiple methods harbored genes for the growth and development of wheat root and shoot organs. For instance, QTN locus Chr2B: AX-94457792 for SRL is comprised of genes TraesCS2B02G406800 and is predicted to encode Formin-like protein which regulates root-hair elongation in rice [32]. This locus also controls the actin cytoskeleton in root hair in wheat [33]. The other largely effected QTN associated with SL (Chr7BS: AX-94876335) was linked to a transcriptional responses of maize seedling root to phosphorus starvation [34]. Most QTNs associated with SRN were at loci containing genes predicted to encode Histone-lysine N-methyltransferase ATXR2, ATPase, F1/V1/A1 complex, alpha/beta subunit and the N-terminal domain superfamily. Other harboring genes encode nucleotide-diphospho-sugar transferases, ribosome-inactivating protein superfamily and papain-like cysteine peptidase superfamily. The large effect QTNs Chr7AS: AX-95249973 and Chr7BS: AX-95123855 associated with FOLRN were involved in plant growth, development and adaptation (Supplementary Table S6).

Discussion
Root traits play a major role in resource uptake, particularly with respect to water access and uptake in the context of drought tolerance. In this way, root traits help maintain crop yield under limited water environments. To realize maximum yield in wheat, a welldeveloped root system is warranted [28,35]. The discovery of a novel source of germplasm and new alleles to improve RSA and introgression of new traits into adapted but otherwise susceptible phenotypes are a desirable approach toward breeding for drought tolerance. Different strategies were adopted for early screening of the wheat RSA, assuming that genotypes with a diverse nature of root architecture at the seedling stage would also respond in a similar way at the adult stage when water and/or nutrients become scarce for grain yield [36]. We studied the root and shoot behaviors at the seedling stage under controlled conditions in perlite vermiculite mixtures in pots. This method of screening was considered a reliable method to examine the root system variations compared to field conditions, where there are several confounding impacts and extracting and measuring roots is difficult. Genetic variability in root systems during seed germination and seedling development is a key trait for seedling establishment and early vigor under water stress and resource poor environments. Improved root traits during this key stage will impact crop emergence and seedling establishment resulting in a good plant population which is one of the key components of yield under both optimal conditions and stress conditions.

Phenotypic Variability
The significant variation observed for various RSA features and shoot traits could be attributed to diversity among the genotypes of the association panel due to their diverse genetic background and wide geographical distributions. The association panel genotypes were collected from various parts of India and around the world. Further, a broad range of variation was detected for individual RSA traits. The highest and lowest morphometric values for various RSA traits (Table 1) Table S2). Li et al. [37] reported moderate to high H 2 for root and shoot traits, ranging from 56.0% to 94.6%.
For acquisition of water and nutrients from the deep zone of soil, RV with a greater number of root tips is essentially required [8,10]. Further, it would be desirable to have genotypes with longer RV, RDW, and other associated parameters so they could extract water and maintain growth under moisture stress environments. The genotype IC406521 from Uttarakhand (India) showed the highest mean value of RV (0.35 cm 3 ), and lowest AD (0.31 mm), thus its root distribution occupied a larger area of soil. Further, the geno-type IC539574 ranked one with the value of 240.66 for FOLRN whereas the genotype IC82425A recorded the highest mean value of 490.33 for SOLRN. Similarly, top ranking genotypes for RLD, LRS, LRD and AD were IC542076 (0.75 cm −2 ), IC82425A (0.86), IC29008 (5.89 cm −1 ) and EC187159 (0.49 mm), respectively. These genotypes could serve as donors of these traits for crop improvement to combat drought stress. The genotypes with the highest AD are considered best for drought stress tolerance due to large xylem vessels which can extract more water and nutrients from deep soil layers [12,38]. We recorded variation for SRN in our diverse panel, ranging from 2 to 6 which was also observed by Djanaguiraman et al. [14]. However, Cane et al. [36] reported 4 to 6 seminal roots in elite wheat varieties which might be due to less variability. We identified donors for multiple traits. For example, the genotype EC426644, an Australian cultivar (Tincurrin), showed high performance for traits, namely TRS, RLD, FOLRN and RSDWR. Bustos et al. [39] also reported Tincurrin cultivar with high root system size in their study.

Multivariate Analyses
Correlations between root and shoot traits shows the balance between the organs of roots and shoots and resource partitioning between the above ground and below ground plant parts [40,41]. In this case, FOLRN and SOLRN exhibited strong correlations with all the traits except AD and SRL. FOLRN directly or indirectly contributes to the variability in root morphological traits which increase water use efficiency at critical stages of plant growth through more absorption of water and nutrients from soil sub layers. The inverse relationship between AD with other traits observed in the present study was in line with the previous studies on flax (Linum usitatissimum L.) and Arabidopsis thaliana [42,43]. This is because fine roots are the main root feature, and the number of root tips increases the root surface area and volume, hence enhancing the absorption efficiency of water and nutrients [8,10]. A high correlation between RV and RDW is obvious. Moreover, a high correlation between shoot and root dry weights observed might be due to supply of nutrients from root to the shoot parts, as is evident in the case of rice [44]. One of the easily scorable key traits, SRN, showed a high correlation with a majority of the root traits, suggesting its utility to provide a broad idea about the root system. The high correlation between root and shoot traits is due to reliance of crops on the root system for water and nutrient uptake.
Based on PCA, RV, RDW, SOLRN, RLD, TRS, RL, LRS, FOLRN, SRN and RSDWR were major contributors to PC1 and most significantly associated with genotypic variations. Likewise, SDW and SL are major contributors to PC2. Hence these traits need to be focused on root studies. RV indicates the proliferation of roots which is essential for plants to uptake more water and nutrients from the soil, whereas linear root elongation may be quite useful in extracting water from deeper soil layers in the event of drought. However, RDW varies among genotypes with varying RL. Therefore, it would be desirable to have genotypes with high RL as well as RDW for a more efficient root system. Higher root density is found to be useful for the plants to uptake more nutrients and therefore the plants with a deeper root system would be more advantageous for drought tolerant genotypes [45,46]. Clustering of genotypes in six groups was also evident in a two-dimensional PCA biplot based on root and shoot morphological traits. Genotypes from two contrasting clusters 1 and 6 could be used in crossing program for genetics studies, mapping population development and trait introgression.

Genome Wide Association and Candidate Genes Identification
Identification of a novel source of genes and genomic regions associated with RSA traits at the seedling stage is expected to boost development of high yielding drought tolerant varieties. In this context, historical wheat germplasm collections maintained in gene banks could prove valuable genetic resource for searching genes/genomic regions using association genetics approach. In this study, we used a diverse panel of 136 wheat genotypes including landraces, varieties, local collection, synthetic germplasm and elite breeding line for mapping RSA and shoot traits. For conducting GWAS, ML-GWAS models were used as these are considered superior to SL-GWAS models for mapping complex traits due to the fact that in ML-GWAS models, all-marker effects are simultaneously estimated. Moreover, unlike SL-GWAS models, these do not require testing of identified associations using stringent multiple testing corrections that generally result in rejection of significant associations [47]. Using six different ML-GWAS models, 456 QTNs were detected for 15 RSA features and shoot traits, distributed over 21 chromosomes. The phenotypic variation estimated (PVE) (r 2 ) for the traits ranged from 0.12% to 38.67%, indicating that RSA traits are controlled by multiple loci with small to moderate effects. This also revealed the complex genetic control of these traits at an early stage of crop growth. Among the six ML-GWAS models used in our analysis, the pLARmEB model was the most powerful, which revealed the maximum number of associations (194) whereas the FASTmrEMMA was the least powerful as it detected the lowest number of MTAs (18). These findings were consistent with the observation of the Safdar et al. [48] who used ML-GWAS models to dissect agronomic traits in bread wheat.
Among the root traits, TRS was found to be associated with drought tolerance in wheat due to spreading of roots in the soil and its effect on the resource uptake [49]. The SNP marker AX-94952472, AX-95105488 and AX-95119337 on chr2AS, 4AL and 7BS were associated with TRS as well as RLD with phenotypic variance in the range of 3.29-11.23%. Root length is an important parameter that determines the ability of a plant to capture water from deeper soil layer [50]. Eight QTNs for RL were detected with more than two models, out of which Q.RL-6BL was consistently detected in five models. Further, this marker was pleiotropic and was associated with RV and RLD which were positively associated traits. Two reliable QTNs for RL (Q.RL-6BL and Q.RL-7AL) were not documented in earlier studies. There, they are considered novel QTNs, whereas other QTNs were reported to be associated with wheat seedling root development under abiotic stresses and phosphorus starvation [51,52]. For traits, FOLRN and SOLRN, two and nine QTNs were identified, respectively. Lateral roots play an important role for foraging water from a shallow depth. Wheat cultivation in the dry area requires more lateral roots at deeper layers for absorption of nutrients and water from deep soil layers. Our analysis revealed a pleiotropic marker, AX-95123855 on 7BS associated with RSA traits FOLRN, SOLRN and LRS. This marker was located within a gene (TraesCS7B02G086500) that encodes Calcineurin which was reported to control inward ion flux in the root, a process essential for plant survival and growth [53].
The SRL is a widely used morphological parameter for root trait which characterizes economic aspects of root system, demonstrating the role of root mass for nutrient acquisition. Genotypes with large SRL have thin roots that increase the surface area per unit root volume enhancing water and nutrient uptake efficiencies [54]. Furthermore, genotypes with small root diameter and large RLD at depth are better adapted to drought conditions [55]. Our study identified 2 QTNs for SRL, one each on 5AL and 2BL. Of these, QTN on 2BL (Q.SRL-2BL) was associated with a gene encoding for Formin-like protein, a family of actin organizing protein involved in root hair development [33]. Further, RLD is considered a primary driver of drought avoidance and enables complete uptake of the soil moisture. Breeding for plants with increased RLD in medium and deep layers and less RLD in shallow soil layers have been proposed as an efficient growth strategy where deep water could be available to crops during late maturity [45,[56][57][58]. Nine QTNs were identified for RLD in our study. Root diameter is another critical parameter determining nutrient uptake and transport. In this study, we identified as many as eleven consistent QTNs for AD and all of them had minor effects, suggesting complex genetic regulation of this trait. Among these, two QTNs Q.AD-1BL (AX-94620468) and Q.AD-6AL (AX-95244609) were detected in three or more models and could represent important genomic regions controlling this trait. Moreover, we found that marker AX-95244609 for QTN Q.AD-6AL was associated with other traits such as SOLRN and LRS, and thus could be a potential target for improvement of more than one RSA traits.
Another important trait is SRN, which emerges first from coleorhiza of seed embryo. The gravitropic response of wheat seminal roots was heritable and suggested to be under control of a single dominant gene [59]. Soriano et al. [60] reported meta-QTNs for SRN in wheat on 2A, 2B, 3B, 4A, 5A, 7A and 7B. In this study, a total of six novel QTNs for SRN were identified with the QTNs on 2D (Q.SRN-2DL: AX-94551988) having consistency in four methods. We found that SNP marker AX-94551988 associated with Q.SRN-2DL was located within gene TraesCS2D02G478700 encoding for histone-lysine N-methyltransferase. In a previous study on Arabidopsis, histone-lysine N-methyltransferase gene was reported to inhibit lateral root development [61]. Based on this, it could be suggested that identified Histone-lysine N-methyltransferase candidate gene might participate in regulating SRN through modification of histone proteins.
Moreover, root biomass accumulation at the seedling stage is beneficial for drought escape as it could enable efficient uptake of water and nutrient from soil. In previous studies, RDW QTLs have been identified under both normal and drought stress conditions [60]. Among the four QTNs for RDW identified in this study, the one on 7AS (Q.RDW-7AS) explained maximum percentage of phenotypic variance (18.35-18.60%) and coincided with a meta-QTL, root_MQTL_68 [60], suggesting this is an important genomic region associated with root trait and its role can be further validated. Among the studied traits, RSDWR is considered the most important as it indicates ability of a genotype to absorb water and nutrients from the soil. Generally, plants have higher root shoot ratio under nutrient deficit soils, suggesting more allocation of resources for the development of root so the plant can extract nutrients from deep and wider zone [62]. Studies have identified genomic regions for the RSDWR trait on various chromosomes of wheat under both control as well as nutrient deficit conditions [60]. Our analysis revealed three QTNs for RSDWR, one each on 2A, 2B and 3B respectively, however each of them have small effect, indicating limited potential for their utilization in breeding programs.
The underground root traits significantly affect development of above ground traits including shoot traits such as, SDW and SL, and thus grain yield. We could identify six and five reliable QTNs for SDW and SL, respectively. The strongest QTN for SDW was identified on 4A, Q.SDW-4AL which explained phenotypic variance in the range of 8.09-19.05%, and was found reliable, as detected by using all five models. The SNP marker associated with this QTN was located within a gene TraesCS4B02G095100 that encodes for an F-box-like domain superfamily. F-box family proteins play a diverse role in plant growth and development and could be critical for shoot growth as their expression has also been detected in the shoot tissues of plant species [49].
Applications of root phenotype-genotype association through GWAS has enabled the identification of important QTNs for root traits that impact shoot traits including yield [63,64]. Markers significantly associated with these traits can be used in marker assisted backcross breeding for varietal improvement. Pleiotropic loci with consistent effects should be amenable to MAS for many traits together. Further, our study provides understanding into phenotype-genotype associations for early root and shoot traits of diverse wheat genotypes by identifying QTNs and proposing plausible candidate genes for future investigations.

Experimental Material and Design
Phenotyping of 140 wheat germplasm of diverse panel including mini core subset was carried out at the Indian Council of Agricultural Research-National Institute of Plant Biotechnology (NIPB), New Delhi (India) as pot screening in growth chambers under controlled environmental conditions during Kharif season (June-September), 2019. Detailed information about the origin and pedigree of the material is given in Supplementary Table S1. Root morphology and root system architecture (RSA) along with shoot traits were studied by destructive methods at early seedling stage (Fifteen-day-old seedlings). Figure 6 describes the steps and methodology used in this experiment. The experimental design adopted for screening was a completely randomized design (CRD) with three replications including two checks (C306 and HD2967). For each replication five biological replicates were used in each pot. The experiment was done in batches for 20 genotypes at each time. Before germination, 20 seeds of each accession were washed carefully and thoroughly with double distilled water and surface sterilized with 0.5% Sodium hypochlorite solution for 30 s. The sterilized seeds were then washed with double distilled water three to four times to remove any trace of adhering chemicals. The seeds were placed well spread in a thoroughly moist germination paper/filter paper taken in a petri dish and allowed to germinate under a growth chamber at 22 ± 1 • C room temperature in the dark.
In this experiment, five germinated wheat seeds were sown initially in pots (4-inch diameter) containing perlite and vermiculite (1:2 ratio v/v). Germinated seeds that were two days old with primary roots of 1 cm in length were placed inside the pot with the root radical facing down at equal depth. The Murashige-skoog (MS) liquid media were supplemented with 8 mM NO 3 − as nutrient media at uniform intervals during the entire growth period. The calcium concentration in all media was kept constant by adjusting the amount of CaNO 3 . Murashige-skoog (MS) N + nutrient solution was prepared as per the given composition and stored in 4 • C. For one liter of working solution preparation, 100 mL of MS media, 4 mL of CaNO 3 and the remaining 896 mL of distilled water constituted 1 L of solution, which were added in equal proportion to all the genotypes grown in controlled conditions at equal intervals. In the growth chamber, 150-200 µmol photon/m 2 /s light intensity, 10/14 dark/light hours, 70% relative humidity (RH) and at 22 ± 1 • C conditions were maintained as described by Sinha et al. [65].

Root and Shoot Traits Measured for Phenotyping
After completing 15 days of growth in the controlled growth conditions, the seedlings were harvested from the pots by removing the perlite + vermiculite mixture and plants were separated carefully without damaging the roots and shoots. The plants were immediately placed on a tray to wash out the adhering perlite + vermiculite particles by running tap water with the utmost care so that to avoid breakage of any roots. The entire seedling was carefully spread on blotting paper and maximum root length and shoot length were measured with the help of a length measuring scale. The root was then extracted intact by cutting at the collar region using a sharp blade and the roots were placed in a tray containing distilled water. The roots were then individually scanned in an Epson Perfection V 700 Photo ® flatbed scanner at a resolution of 400 dpi modified for this purpose (Regent Instruments Inc., Quebec, QC, Canada) as per the manufacturer's guidelines. The root images from the scanner were analyzed with customized software WinRHIZO™ (Regent Instruments Inc., Quebec, QC, Canada) [66]. Various root trait data were recorded by software which were later transformed in major RSA traits measured manually using published protocols [32,[67][68][69] (Table 1).

Statistical Analyses of Phenotypic Data
Mean data across the replications were used as the input data for statistical analyses and GWAS analysis. Descriptive statistics and frequency distribution were analyzed to check range of variability among the traits. Pearson's correlation coefficient, cluster analysis and PCA was performed using SAS software version 9.3 (JMP) program (SAS Institute, Cary, NC, USA). Heritability (H 2 ) was estimated from the analysis of variance.

DNA Extraction and SNP Genotyping
Genomic DNA was extracted from 15-day-old wheat seedlings using CTAB methods [70]. DNA quality was checked using Nanodrop TM 2000 (Thermo Fisher Scientific, Wilmington, DE, USA). Samples with good quality DNA were genotyped using Axiom ® Wheat Breeders' Array (Thermo Fisher Scientific, Wilmington, DE, USA) according to the procedure described by Affymetrix (Axiom ® 2.0 Assay for 384 samples P/N 703,154 Rev. 2). SNP markers with >10% missing data and <10% MAF were excluded. Identified SNPs were localized on a wheat genome assembly International Wheat Genome Sequencing Consortium (IWGSC) RefSeq version 1.0 using BLASTn program with default parameters.

Population Structure and LD
Population structure was estimated with 525 unlinked SNP markers nearly uniformly distributed across the wheat genome using the Bayesian model-based approach implemented in STRUCTURE program version 2.2 (Pritchard Lab, Stanford University, Stanford, CA, USA) [71]. A burn-in of 20,000 iterations followed by 50,000 Monte Carlo Markov Chain (MCMC) was run to estimate the number of subpopulation (k) in a putative range of k = 1 to 10. The subpopulation number was estimated using an ad hoc statistic delta k based on the rate of change in log probability of data between successive values [72]. The squared allele frequency correlation (r 2 ) between SNP markers was used to estimate LD across sub-genome A, B, and D using TASSEL v5.0 (Buckler's Lab, Ithaca, NY, USA) [73]. LD decay across the three sub-genomes and whole genome level was estimated as the physical distance between SNPs where average r 2 reduced to half of the maximum LD value.

Genome Wide Association Analysis
The use of multi-locus methods that capture small effect loci in complex polygenic traits such as in plant roots and shoots have recently become a popular and feasible approach. To benefit the algorithmic merits of different models and support results of one by another, it is also advantageous to apply multiple methods [49]. GWAS analysis was performed using six multi locus GWAS methods within mrMLM [74], FASTmrMLM [75], FASTmrEMMA [76], pKWmEB [77], pLARmEB [78] and ISIS EMBLASSO [79], which were included in the R package mrMLM v3.1 [80]. All parameters were set at default values in this GWAS. The critical thresholds of significant association for the six methods were set as LOD score 3.00 or >3.00. The most significant QTNs, detected in at least two methods, were considered as reliable QTNs. The associated SNPs and their putative underlying genes were illustrated on the wheat chromosomes using Map Chart 2.3, (https: //www.wur.nl/en/show/Mapchart.htm accessed on 8 June 2021) [81]. The favorability of alleles at QTNs detected by at least four of the multi-locus models was illustrated using a box plot based on mean phenotypic value of genotypes with each allele.

Identification of Potential Candidate Genes
SNPs (probe sequences) that were significantly associated with root architecture traits were searched against the Triticum aestivum genome assembly IWGSC-refseq version1.0 in online web resource Ensemble plants (https://plants.ensembl.org/Triticum_aestivum/ Tools/Blast, accessed on 14 November 2020) using BLASTn with default parameters to identify potential candidate genes. BLAST2GO tool was used to get annotation of expressed transcripts [82].

Conclusions
This study revealed wide variability for RSA and shoot traits at the seedling stage in the studied association panel. Ten top and bottom performing lines for 15 traits were identified for use in genetics, development of mapping population and introgression studies. Based on PCA, several traits such as RV, RDW, SOLRN, RLD, TRS, RL, LRS, FOLRN, SRN and RSDWR, were the most influential traits for phenotypic variations. The GWAS analysis enabled genetic detection of RSA traits and revealed 11 pleotropic loci associated with correlated traits. Further, putative candidate genes were identified from the associated genomic region that could be validated using a functional genomics approach. Development of wheat cultivars possessing superior root traits will play an important role in enhancing drought tolerance under water stress conditions. The large variability found in Indian germplasm for RSA traits and the novel genomic regions regulating them makes this germplasm a valuable source for improving root architecture, which plays a significant role in absorption and uptake of water and nutrients and increases crop productivity.