Genetic Dissection Uncovers Genome-Wide Marker-Trait Associations for Plant Growth, Yield, and Yield-Related Traits Under Varying Nitrogen Levels in Nested Synthetic Wheat Introgression Libraries

Nitrogen is one of the most important macronutrients for crop growth and metabolism. To identify marker-trait associations for complex nitrogen use efficiency (NUE)-related agronomic traits, field experiments were conducted on nested synthetic wheat introgression libraries at three nitrogen input levels across two seasons. The introgression libraries were genotyped using the 35K Axiom® Wheat Breeder's Array and genetic diversity and population structure were examined. Significant phenotypic variation was observed across genotypes, treatments, and their interactions across seasons for all the 22 traits measured. Significant positive correlations were observed among grain yield and yield-attributing traits and root traits. Across seasons, a total of 233 marker-trait associations (MTAs) associated with fifteen traits of interest at different levels of nitrogen (N0, N60, and N120) were detected using 9,474 genome-wide single nucleotide polymorphism (SNP) markers. Of these, 45 MTAs for 10 traits in the N0 treatment, 100 MTAs for 11 traits in the N60 treatment, and 88 MTAs for 11 traits in the N120 treatment were detected. We identified putative candidate genes underlying the significant MTAs which were associated directly or indirectly with various biological processes, cellular component organization, and molecular functions involving improved plant growth and grain yield. In addition, the top 10 lines based on N response and grain yield across seasons and treatments were identified. The identification and introgression of superior alleles/donors improving the NUE while maintaining grain yield may open new avenues in designing next generation nitrogen-efficient high-yielding wheat varieties.

Nitrogen is one of the most important macronutrients for crop growth and metabolism. To identify marker-trait associations for complex nitrogen use efficiency (NUE)-related agronomic traits, field experiments were conducted on nested synthetic wheat introgression libraries at three nitrogen input levels across two seasons. The introgression libraries were genotyped using the 35K Axiom ® Wheat Breeder's Array and genetic diversity and population structure were examined. Significant phenotypic variation was observed across genotypes, treatments, and their interactions across seasons for all the 22 traits measured. Significant positive correlations were observed among grain yield and yield-attributing traits and root traits. Across seasons, a total of 233 marker-trait associations (MTAs) associated with fifteen traits of interest at different levels of nitrogen (N0, N60, and N120) were detected using 9,474 genome-wide single nucleotide polymorphism (SNP) markers. Of these, 45 MTAs for 10 traits in the N0 treatment, 100 MTAs for 11 traits in the N60 treatment, and 88 MTAs for 11 traits in the N120 treatment were detected. We identified putative candidate genes underlying the significant MTAs which were associated directly or indirectly with various biological processes, cellular component organization, and molecular functions involving improved plant growth and grain yield. In addition, the top 10 lines based on N response and grain yield across seasons and treatments were identified. The identification and introgression of superior alleles/donors improving the NUE while maintaining grain yield may open new avenues in designing next generation nitrogen-efficient high-yielding wheat varieties.

INTRODUCTION
The global demand for nitrogen currently stands at about 117 million metric tons with a projected annual increase of approximately 1.5% expected in the near future (FAO, 2019). Farmers generally apply high doses of nitrogenous fertilizers to ensure good yields. The high input of commercially available fertilizers has led to the degradation of air, soil, and water quality (Hickman et al., 2014;Russo et al., 2017). In addition, when the supply of nitrogen (N) is in excess of crop N demand, it increases the susceptibility of plants to various diseases and insect pests (Reddy, 2017). Therefore, it is necessary to optimize and improve the nitrogen use efficiency (NUE) of cereal crops to maximize yield in addition to minimizing the negative impact of increase in N use on the environments and natural resources. Identification of marker-trait associations (MTAs) can be used to make effective targeted introgressions and is one possible genetic method to address the challenge of developing N-efficient wheat varieties with stable yield under N-limited environments.
Wheat varieties that maintain yield under moderate or intense N deficiency can adapt to low input systems. To breed such varieties, genetic variation for adaptation traits to N deficiency is required. To date, limited quantitative trait loci (QTL) for both yield and its response to N deficiency in wheat under field conditions have been documented. Detection of genotypes and underlying QTLs for maintaining yields at low N levels are of value in wheat breeding programs designed to increase N-deficiency tolerance. Some QTLs influencing N uptake have been genetically mapped in wheat under different doses of fertilizer application using biparental populations Laperche et al., 2008;Xu et al., 2014;Mahjourimajd et al., 2016;Deng et al., 2017). A number of genetic loci for agronomic traits related to N use and grain yield have also been mapped to the chromosomal regions containing the GS2 gene in wheat and rice (Prasad et al., 1999;Yamaya et al., 2002;Obara et al., 2004;Habash et al., 2007;Laperche et al., 2008;Fontaine et al., 2009). This suggests the role of the genomic region surrounding GS2 (Pritchard and Wen, 2004) is favorable in breeding wheat and rice varieties with improved agronomic performance and nitrogen use efficiency (NUE). Other genetic regions associated with N uptake have also been detected in rice (Wissuwa et al., 1998;Ming et al., 2000), wheat (Su et al., , 2009; maize (Zhu et al., 2005), common bean Yan et al., 2004), and soybean (Li et al., 2005;Liang et al., 2010). The NRT2.1, NRT2.2, and NAR2.1 genes have been reported to be the important contributors to the high-affinity transport system in Arabidopsis roots (Orsel et al., 2006). Sixteen genes were identified in wheat homologous to characterize the low-affinity nitrate transporter NPF genes in Arabidopsis, suggesting a complex wheat NPF gene family (Buchner and Hawkesford, 2014). The regulation of wheat NPF genes by plant N-status indicated the involvement of these transporters in the substrate transport in relation to Nmetabolism.
The phenotypic traits reported to be associated with NUE in cereal crops so far include root number, length, density, and branching (Morita et al., 1988;Yang et al., 2012;Steffens and Rasmussen, 2016), dense and erect panicle (Sun et al., 2014), plant height (Gaju et al., 2011), and leaf width (Zhu et al., 2020). The collocation of QTLs for N-uptake and root architecture traits have suggested that breeding for better and efficient root systems is a way to improve NUE (Coque et al., 2008;Sandhu et al., 2015).
Diverse accessions, landraces, breeding populations, and nextgeneration mapping populations, including nested-association mapping (NAM) and multi-parent advanced generation intercross (MAGIC) populations have shown a potential for mining novel genetic variation in rice Sandhu et al., 2019;Subedi et al., 2019), wheat (Mackay et al., 2014), maize (Yu et al., 2008), and soybean (Xavier et al., 2015). The NAM and MAGIC populations have proven advantageous over biparental populations as they capture additional recombination breakpoints, thus increasing the allelic diversity and improving the power of QTL detection (Yu et al., 2008;Scott et al., 2020). Further, the availability of high throughput genotyping platforms to generate uniformly distributed genome-wide molecular markers is critical for the high-resolution genetic dissection of polygenic traits, and the tracking of favorable alleles in breeding populations (Pandey et al., 2012(Pandey et al., , 2016Varshney et al., 2013). To date, a series of high-density wheat single nucleotide polymorphism (SNP) arrays, such as the Illumina 9K iSelect SNP array (Cavanagh et al., 2013), Illumina 90K iSelect SNP genotyping array (Wang et al., 2014a), 15K SNP array (Boeven et al., 2016), Axiom R 660K SNP array, 55K SNP array, Axiom R HD 820K genotyping array (Winfield et al., 2016), 35K Axiom array (Allen et al., 2017), and 50K Triticum TraitBreed array (Rasheed and Xia, 2019) have been developed and their utility has been demonstrated across a range of applications.
In the present study, we developed nested synthetic wheat introgression libraries capturing novel genetic variation. The libraries were genotyped using a high-density SNP array and phenotypically assessed for root traits and agronomic performance under three N input conditions in the field. Genome-wide association mapping was used to identify MTAs for the root and agronomic traits, and lines carrying favorable genetic combinations were also identified for use in future breeding for improved N use.

Agronomic Practices and Management of Experiments
The N-SHW library, six parents, and two synthetic hexaploid wheats were assessed under field conditions at the experimental farms of School of Agricultural Biotechnology, PAU Ludhiana (30 54' N latitude, 75 48' E longitude, and 247 m above sea level) over 2 years in 3N levels (6-year x N combinations). The soil analysis of experimental plots showed soil pH: 7.3; EC: 0.157 ds/m; organic content: 0.29%; N: 99 Kg ha −1 , P: 69 Kg ha −1 ; K: 116 Kg ha −1 , Zn: 4.77 parts per million (ppm), Cu: 1.73 ppm, and Fe: 7.30 ppm. The N-SHW libraries were evaluated in rabi seasons of 2018 and 2019. Details of the number of lines tested and experimental design is provided in Table 1. The breeding material was sown on 21 st of November and 18 th of November in 2018 and 2019, respectively. In both the years, the experiments were conducted at three N levels [i.e. zero N (0 Kg ha −1 ), half N (60 Kg ha −1 ), and full N (recommended, 120 Kg ha −1 ),] referred to as N0, N1, and N2, respectively. The recommended dose of phosphorus, potassium, and manganese was applied at the time of sowing. Half of N was applied at the time of sowing while the other half was applied in two equal splits, the first at crown root initiation stage and the remaining at the maximum tillering stage in both the N1 and N2 experiments. N0 was treated as a control. Recommended fungicides and insecticides were applied to control stripe rusts, brown rusts, and aphids at jointing, booting, and 10 days after anthesis to prevent diseases and pests. Weeds were controlled manually.

Characterization of Phenotypic Traits
Under field conditions, a total of 22 traits were assessed in all experiments across both seasons except the maximum root length and root angle which were measured in 2018 only. The details of  Shoots were separated from the roots, and the fresh root weight (FRW; g) and fresh shoot weight (FSW; g) were measured. The root and shoot samples were dried at 70 • C in an oven until constant dry root weight (DRW; g) and dry shoot weight (DSW; g) were observed, while the roots were cleaned thoroughly and stored in 70% alcohol at 4 • C for root trait evaluation. The maximum root length (MRL) and root angle (RA) were measured using ImageJ software. total root length (TRL), total root surface area (RSA), total root diameter (RD), total root volume (RV), number of forks (NF), and number of tips (Ntips) were recorded using WinRhizo STD4800 (Supplementary Figure 2). The roots were then dried at 70 • C in the oven until constant DRW was observed. The data on N-uptake related traits were recorded using a chlorophyll meter, Soil Plant Analysis Development Meter (SPAD502) and leaf color chart (LCC). The LCC provides a decision-support system to the farmers for sustaining the high yields with optimum N dose in the field crops. It measures the leaf color variations of 6 SPAD units comprising 3, 3.5, 4.0, 4.5, 5.0, and 6.0 and provides N recommendation in the field crops. Flag leaf length (FLL) and flag leaf width (FLW) were recorded using a centimeter scale. Days to flowering (DTF) rate of about 50% was recorded when 50% of the plants in a plot exerted their panicles. Spikelets per spike (SPS) were counted manually from five random plants. The number of productive tillers (NPT) was counted manually in 0.5 m row length and shoot biomass (SB) at harvesting was measured from 0.5 m row length. The plant height (PHT) in cm was measured as the mean height of five random plants for each entry measured from the base of the plant to the tip of the panicle at the maturity stage. The plants were harvested at physiological maturity or when 80-85% of the panicles turned to golden yellow and the panicles at the base were already at the hard dough stage; harvested grains were threshed, dried, and weighed to determine the grain yield (GY). The shoot biomass (SB) was recorded at harvesting.

Phenotypic Data Analysis
Analysis of variance, and experiment-wise mean for each season was calculated using a mixed model analysis in PBTools V 1.4.0. for augmented design and in STAR Version: 2.0.1 for the split plot design. In the split plot design, the N levels were considered as the main plot and the breeding lines as subplot. Fisher's ttest was used to determine the significant difference among the breeding lines, treatments, and to estimate the interactions. The correlation analysis among traits was performed in R. v.1.1.423.
To evaluate the phenotypic stability and GY adaptability of the breeding lines across seasons and treatments, the genotype and genotype × environment (GGE) biplot analysis was performed, considering the effects of genotype (G) and genotype by environment (GE) as random. The best linear unbiased prediction (BLUP) values of the G and GE effects were calculated. The multiplicative model in PB tool version 1.3 (http://bbi.irri.org/) was used to explain the relationship between G and the seasons.

Genotypic Data
High-density genotyping was performed using the 35K Axiom R Wheat Breeder's Array (Affymetrix UK Ltd., United Kingdom). The quality pre-processing of 35,143 markers obtained from the 35K chip was done using PLINK software (Purcell et al., 2007). A total of 9,474 single nucleotide polymorphisms (SNPs) with minor allele frequency (MAF) of >5%, maximum heterozygote proportion of 0.1 and missing rates <0.1 were used to estimate the genetic relationships and for the mapping of MTAs for different traits associated with plant growth, yield, and yieldrelated traits. Principal component analysis (PCA) was carried out to detect and correct for population structure.

Population Structure and Association Analysis
The model-based STRUCTURE V. 2.3.4 software was used to test K values from 1 to 10, with a burn-in period of 10,000 and 1,000,000 Markov chain Monte Carlo (MCMC) reps after burnin in order to assess the population structure in the 352 breeding lines using a total of 9,474 SNPs. The consistency and accuracy of the results were validated across 10 runs for each K. The K value with maximum likelihood over the 10 runs was used to estimate the most appropriate number of clusters (Pritchard and Wen, 2004). The population structure was determined by plotting the proposed number of subpopulations against the delta k (Earl, 2012). The PCA was performed in Genome Association and Prediction Integrated Tool (R/GAPIT) and added iteratively to the fixed model, ranging from PC1 to PC10.
Significant marker-trait associations were identified using a compressed mixed linear model (CMLM)/P3D (population parameters previously defined) in GAPIT executed in R. Identity by state (IBS) values and a relatedness matrix were used to estimate the random effect and genetic similarity of the accessions, respectively. The statistical power of the association studies was further improved by considering the population structure (Q value) and kinship matrix (K) estimated from the genotyping data. The Bonferroni correction method was used to correct for false positives in the analysis, using the stringent pvalue benchmark. The Bonferroni multiple test correction was performed (0.05/9474; significance level of 5%/total number of markers used in analysis) and the calculated threshold value was found to be 5.28 × 10 −6 . The allelic effect of all the significant markers associated with the measured traits was determined by comparing the mean phenotypic values and the significant allelic variants for the trait/s using a Kruskal-Wallis test in R.

Candidate Gene Analysis and Functional Annotation of Putative Candidate Genes
Single nucleotide polymorphisms that exhibited a false discovery rate (FDR) corrected p < 0.05 for a particular trait of interest were evaluated as markers for the potential putative candidate genes. A window of 1 Mb adjacent to each significant SNP was examined for candidate genes and annotations were identified through the Ensembplants database (http://plants.ensembl.org/index.html).
The functional annotation and gene ontology (GO) of the identified putative candidate genes were performed using OMIX box software. Blast (E ≤ 1 0 −5 ) was performed using the CloudBlast tool against Triticum (nr_subset) [monocots_triticum, taxa:4564] and NCBI non-redundant database (http://www.ncbi.nlm.nih.gov), followed by the InterPro using CloudIPS, followed by GO mapping, and annotation configuration. The GO terms were then used to generate the semantic similarity-based scatterplots/interactive graphs/tag clouds by using REVIGO (http://revigo.irb.hr/).

Defining N-Irresponsive and N-Responsive Lines
The genotypes that showed more or equal/stable GY with the minimal application of nitrogen fertilizer when compared to the recommended or standard N fertilizer application were considered as the N irresponsive (NIR) genotypes or the top grain yielders across seasons and treatments. On the other way around, the genotypes that were low yielding or not able to maintain the GY with the minimal application of N fertilizer when compared to the recommended or standard N fertilizer application, were considered as the N responsive (NR) genotypes or the poor grain yielders across seasons and treatments.

Significant Phenotypic Trait Variation and Correlations Detected Across Nitrogen Treatments
The 352 N-SHW lines, six parents, and two synthetic hexaploid wheat donors were screened for 22 traits in six growing conditions (2 years x 3 nitrogen levels). The ANOVA revealed significant genetic variation for the root, plant morphological, and agronomic traits among genotypes, treatments, seasons, and their interactions (genotype x treatment, genotype x season, treatment x season and genotype x treatment x season) ( Table 2). The detailed information on trial means, least significance difference (LSD), and heritability for all the traits measured are presented in Supplementary Table 1 Table 2). The phenotypic data of the traits measured in the present study were averaged across two seasons and are presented as mean values in Supplementary Table 3.
Grain yield increased with applied N level. In the N0 treatment, the average GY of the tested breeding lines across seasons was 2,022 kg ha −1 and it ranged from 564 to 4,092 kg ha −1 (Supplementary Table S3). In the N60 treatment, the GY varied from 882 to 4,685 kg ha −1 with an average GY of 2,357 kg ha −1 and, while in N120 treatment, the GY varied from 1,332 to 4,270 with an average of 2,579 kg ha −1 (Supplementary Table 3). Across seasons, N in the limited conditions (N0) resulted in the 14% and 22% GY reduction compared to N60 and N120 treatments, respectively. The N application also significantly increased the shoot biomass (SB) by 8% in N60 and 52% in N120 treatment across seasons. The average NPT across experiments was higher in N120 (28) compared to N60 (24) and N0 (22) (Supplementary Table 3). Under the N0 treatment, the average value of the leaf color chart (LCC) varied from 3.3 to 4.8, from 3.5 to 5.1 in the N60 treatment, and from 4.1 to 5.3 in the N120 treatment (Supplementary Table 3). The response of lines in terms of average dry shoot weight (DSW) across seasons increased from 3.28 in N0 to 3.62 in N60 to 3.75 in N120 treatment (Supplementary Table 3). The minimum and maximum values of dry root weight (DRW) under N0 were 0.187 g and 2.425 g; 0.298 g and 2.001 g under N60, and 0.338 g and 2.333 g under N120, respectively (Supplementary Table 3). The average root diameter was the highest under N60 (0.610 g) compared to N0 (0.560 g) and N120 (0.409 g) (Supplementary Table 3). Across seasons, the average flowering was delayed by 2 days under the N0 treatment compared to the N60 and N120 treatments. Average plant height (PHT) was lower (92 cm) in N0 compared to N60 (95 cm) and N120 (99 cm).
We calculated the Pearson's correlation coefficients between all the traits measured in N0 (Figure 2A), N60 (Figure 2B), and N120 (Figure 2C) treatments. The Pearson's correlation coefficients across all treatments considering pooled mean data for all the traits measured in the present study is presented in Supplementary Figure 3. The strongest and most TABLE 2 | Analysis of variance (ANOVA) for the NUE-related, root, plant morphological, yield, and yield-related traits among genotype (G), treatments (T), seasons (S), and their interactions (G x T, genotype x treatment; G x S, genotype x season; T x S, treatment x season; and G x T x S, genotype x treatment x season).

Population Structure Analysis Detected Three Genetic Subpopulations
The population structure of the N-SHW lines was assessed to understand the genetic structure of the 352 lines based on 9,474 single nucleotide polymorphisms (SNPs) distributed across all 21 wheat chromosomes. The most appropriate K explaining the population structure was K = 3 at minor allele frequency (MAF) ≥ 5% ( Figure 3A). The kinship heatmap indicated a weak relatedness in the panel ( Figure 3B). The first three principal components (PCs) were the most informative and gradually decreasing (Figures 3C,D) until the tenth PC. The kinship and PCs were considered during the genome-wide association study (GWAS) to correct for population structure. The appropriate number of subpopulations was determined from the largest delta K value of 3 ( Figure 3E).

Mapping Reveals Significant Marker-Trait Associations for All Traits
Genome-wide association study was performed exploiting the phenotypic variability in the 352 N-SHW lines using 9,474 SNPs from the 35K Axiom R Wheat Breeder's Array. Using thelog(P) ≥0.001 at 5% significance level, a total of 233 marker-trait associations (MTAs) were detected across seasons associated with fifteen traits of interest at different N levels (N0, N60, and N120; Table 3). Of these, 45 MTAs were associated with 10 traits in the N0 treatment, 100 MTAs were associated with 11 traits in the N60 treatment, and 88 MTAs were associated with 11 traits in the N120 treatment (Table 3). Across seasons and N treatments, a total of 53 MTAs associated with more than one trait/treatment were detected (Table 3). In addition to these 53 MTAs, another 41 MTAs associated with a single trait only were detected across seasons (Supplementary Table 4). All MTAs detected in the present study either in one season, both seasons, each treatment, or across treatments are compiled in Supplementary Table 5 The Manhattan plots depicting the significant -log (p-values) for the MTAs associated with nitrogen use efficiency (NUE)related traits, root traits, and yield/yield related traits measured in the present study at different levels of N are presented in Supplementary Figures 4, 5 and Figure 4, respectively. Location of significant MTAs and SNP marker density distributed across 21 wheat chromosomes is presented in Figure 5. The SNPs for positively correlated traits, such as GY, BY, tips, RSA, RV, and forks appeared to be collocated on Chr 6A at the different   (Burkhead et al., 2003); signal transduction (Arnesano et al., 2003); nitrogen regulatory response in bacterial and eukaryotic chloroplast (Ninfa and Atkinsonm, 2000;Arcondéguy et al., 2001) AX  Regulates Pi uptake by modulating PHT1;1 expression in Arabidopsis (Wang et al., 2014b); age-triggered leaf senescence (Chen et al., 2017); Benzothiadiazole-inducible blast resistance (Shimono et al., 2007); resistance against F. graminearum in wheat (Bahrini et al., 2011); broad-spectrum resistance to wheat powdery mildew (Cao et al., 2011) AX-95190381 2D 591027900 N0 (FSW,DSW); N60 (FRW, DRW); N120 (DRW,FSW) 2.11E-07 0.18 0 TraesCS2D02G493700 591031752 591027189 Serine-threonine protein kinase Regulates stress-responsive gene expression in rice (Diédhiou et al., 2008), Negative regulator of immune responses in Arabidopsis ; confers durable and broad-spectrum resistance to wheat powdery mildew (Cao et al., 2011) AX-95018936 2D 595159320 N60 (DRW); N120 (DRW,DSW)  Hydrolase activity and root hair cell differentiation (https://www.uniprot.org/ uniprot/Q84JS5); response to salt stress (Jung et al., 2015); involved in male sterility (Mok and Mok, 2001) AX-94962360 2D 596914793 N0 (DRW,FSW,DSW); N60 (FRW,DRW,DSW); N120 (FSW) 4.51E-07 0.185 0.002 TraesCS2D02G503000 596917719 596911132 Plasma membrane H + -ATPase Plant adaptation to environmental stresses (Janicka-Russak and Kabała, 2015), P deficiency and Al toxicity (Wang et al., 2014a;Yu et al., 2015), transport of various nutrients (nitrate, phosphate and potassium) through roots, elongation of hypocotyls in Arabidopsis (Takahashi et al., 2012); NH + 4 metabolism in rice roots (Weng et al., 2020); auxin-mediated cell elongation during wheat embryo development ( Protein transport, response to salt stress (https://www.uniprot.org/uniprot/Q9LFP1), tolerance to water stress (Singh et al., 2017); growth and immune response in Arabidopsis (Yun et al., 2013) AX-94565231 6D 683646420 N0 (RSA); N60 (RSA, TRL,RV,Tips); N120 (RSA,RV,Forks,Tips)      Figure 5). A cluster of 17 SNPs spanning a 7.7 Mb region on the short arm of 6A showed association with GY at N60 and N120 (Figure 5). Across seasons and treatments, significant association in a 198 kb region on the long arm of Chr 6A were detected for root traits (RSA, RV, tips and forks). The SNP, AX-94565231 at 683.64 Mb on the long arm of Chr 6D showed association with different root traits (RSA, RV, tips, and forks) across seasons and treatments. In the N60 treatment, significant associations for FLW were detected in a 5.6 Mb region (549,799,201,748 bp) on the long arm of 4A. Interestingly, the association of the trait, DTF with SNP AX-95136655 on Chr 3B at 234.49 Mb was common under N0, N60, and N120 treatments (Figure 5). In the N0 treatment, a significant association harboring three strongly associated SNPs (AX-94593608, AX-94786978, and AX-95134564) spanning the genomic region 76 bp on the long arm of 3A was detected for NPT (Table 3). Further, single SNPs were identified in association with different traits at different N levels. For example, the SNP AX-94914391 (36.43 Mb, 6A) was significantly associated with SB at N0 and with GY at both N60 and N120 (Figure 5). The SNP 2B) showed association with FSW at N60 and with DRW at N120. The SNP network indicating the significant marker-trait interactions is presented in Figure 6.

Candidate Gene Identification and Functional Annotation
In order to identify candidate genes underlying the consistent MTAs, we surveyed putative candidates in a 1 Mb upstream and 1 Mb downstream region the identified significant SNPs using EnsemblPlants (http://plants.ensembl.org/index.html). Detailed information on the identified candidate genes is presented in Table 3.
The GO term of identified putative candidate genes was categorized into four groups according to their trait relatedness; The NUE uptake-related (LCC, SPAD, FSW, DSW), root morphological [maximum root length (MRL), TRL), RSA, root diameter (RD), RV, NF, N tips, FRW, and DRW], plant morphological [flag leaf length (FLL), FLW, plant height (PHT)] and GY/yield-attributing traits (days to flowering (DTF), SPS, NPT, SB, and GY]. Most of the putative candidate genes in NUE-uptake related traits across treatments were associated with protein phosphorylation/proteolysis, recognition of pollen, molybdoprotein cofactor biosynthetic process, and transmembrane transport (Supplementary Table 6). Some were part of the cellular component organization and molecular functions of binding molecules and ions, catalytic activity, peptidase activity, and transmembrane transport activity (Supplementary Table 6 and Supplementary Figure 6). The putative candidate genes for the root morphological traits were associated with N compound metabolic processes, phosphorylation, proteolysis, catabolic processes, response to stresses, and regulation of flower development by delineating the composition, and architecture of gene regulatory network underlying flower development and carbohydrate metabolism (Supplementary Table 6 and Supplementary Figure 7). The cellular components include chloroplast, ribosome, membrane, cytoplasm, nucleus, and mitochondria (Supplementary Tables 6, 7). The primary molecular functions related to these genes were catalytic activity (protease, peptidase, hydrolase, transferase, ligase, and oxidoreductase), and binding activity (small molecule binding, ion binding, lipid binding, and carbohydrate derivative binding) (Supplementary Table 6 and Supplementary Figure 7). The putative candidate genes for the plant morphological traits were mainly associated with phosphorylation, response to light-intensity, stress-related responses, and metabolic processes. They were related to the molecular functions of metal ion binding, catalytic activity, kinase activity, and DNA/RNA/ATP binding (Supplementary Table 6 and Supplementary Figure 8). The yield and yield-attributing traits-related putative candidate genes were associated with phosphorylation, metabolic process, protein folding, catabolic process, response to water-stress and light, flower development, and pollen recognition (Supplementary Table 6 and Supplementary Figure 8). The molecular functions include catalytic activity (peptidase, hydrolase, lyase, oxidoreductase, and transferase), binding activity (ion, metal, ATP/GTP, polysaccharide, protein, and DNA), and metabolic activity (Supplementary Table 6 and Supplementary Figure 9).

Selection of Promising Lines With Stable Performance for Use in Breeding
To identify stable breeding lines across treatments and seasons, a GGE biplot method was used. The first two principal components (PCs) explained 77.7% (PC1 = 50.3%, PC2 = 27.4%) of the total GGE variation in the data (Figure 7). The ranking of breeding lines based on their mean GY and stability across seasons and treatments (Supplementary Table 7) was used to identify 20 breeding lines with high and stable yields across seasons and treatments (Supplementary Table 8). The percentage increase in the GY of top 20 breeding lines derived from the nested introgression libraries possessing high and stable GY (kg ha −1 ) compared to the respective recipient parent averaged across two seasons under three different N treatments is presented in Figure 8. Based on GY data across seasons and treatments, the top 10 N-irresponsive (NIR-top grain yielders) and 10 Nresponsive (NR-poor grain yielders) breeding lines were also identified (Table 4 and Figure 9).
Further analysis was undertaken to assess the significant differences between the mean values of the allelic classes of MTAs for root growth and GY using the Kruskal-Wallis test. The presence of favorable alleles with significant differences was checked in promising breeding lines. This allowed for the selection of 20 promising breeding lines possessing favorable allele combinations for improving the plant root growth ( Figure 10A) and grain GY under N limitation ( Figure 10B).

DISCUSSION
An increase in crop production by the development of highyielding varieties is largely dependent on the supply of N fertilizers. Excessive application of nitrogenous fertilizer is becoming very expensive which accounts for the great loss of economic profit to the farmers in addition to the negative impacts on the environment (Hawkesford and Griffiths, 2019). The reliable phenotyping under low nitrogen input is very challenging and affected by genotype (G), environment (E), and the G x E interactions (Rao et al., 2018). Proper understanding of the genotype behavior, identification, and development of Nefficient genotypes without compromising the GY is a paramount need for improving the NUE. Notably, very few wheat breeding programs are targeting the development of N-efficient genotypes. In crop plants, such as wheat, the efforts are constrained due to a lack of variation in the cultivated germplasm for NUE. The narrow genetic diversity and fewer recombination events in the biparental mapping populations may result in poor quantitative trait loci (QTL) detection power (Gangurde et al., 2019). The next-generation high-resolution mapping populations, such as nested synthetic wheat introgression libraries used in the present study may provide a vast and untapped source of genetic variations for the NUE-related traits due to high numbers of recombination events. The use of synthetic hexaploid wheat in the present study presents an effective genetic resource for transferring the agronomically important genes from wild relatives to the common wheat. The introgression of favorable alleles associated with root traits and GY from Ae. tauschii wild accessions to cultivated wheat (Figure 10) indicated the potential of synthetic wheat providing new sources for improving yield potential and NUE when bred with the modern wheat varieties.
The different traits associated with N uptake and NUE were studied in nested synthetic wheat introgression libraries at three different nitrogen levels. The ANOVA results revealed the native variation across the genotypes toward the N response which had given the possibility to identify the NUE lines under different levels of N. The genotypic variations purely reveal the phenotypic plasticity of the breeding lines toward traits. The diverse responses have been observed among the breeding lines across different levels of N, despite similar growth conditions and an equal amount of nitrogenous fertilizer application in a given N level as indicated by significant differences among the genotypes within and across treatments and non-significant differences among the replications. Significant G x E, G x S, G x E, and G x T x S interactions indicated that the seasons and environments under different levels of N application were FIGURE 7 | The GGE biplot showing the performance of 352 nested synthetic wheat introgression lines across seasons and treatments (N0, N60, and N120). The environment view refers to the three different levels of nitrogen (N) application: N0, N60, and N120. The genotype view refers to the 352 nested synthetic wheat introgression lines. The numeric number refers to the coding for the introgression lines, which is given in detail in the Supplementary Table 7. a critical factor in explaining the genotypic variance for the traits measured in the present study. The results reported in the present study concurred with other reported studies on rice (Srikanth et al., 2016) and wheat (Sial et al., 2005;Belete et al., 2018).
In general, the increase in GY was correlated with the increase in the rate of N fertilizer application, which might be due to the availability of sufficient N for the proper growth and development of the plants. Šarčević et al. (2014) reported a 10% reduction of GY at low N condition compared to normal condition in wheat. The significant and positive correlation among different root traits and GY and yield-attributing traits indicated complementary functional roles of the root traits in improving the GY by improving the nutrient acquisition from the soil. The collocation of MTAs for the correlated traits strengthens the significance of MTAs. A significant positive correlation between GY and NUE-related traits in wheat, maize, and oilseed rape (Fageria et al., 2010;He et al., 2017;Belete et al., 2018) signified the importance of NUE-related traits in improving the GY under limited N conditions. Different mapping approaches using nested-association mapping (NAM) populations successfully exploited the genetics of complex traits and facilitated the discovery of candidate 4 | The top 10 nitrogen irresponsive (NIR) and 10 nitrogen responsive (NR) breeding lines with contrasting grain yield (GY; kg ha −1 ) derived from pooled mean over two seasons and three treatments. association study (GWAS) keeping into account, the genetic effects produced in each genetic background. The associated SNPs were used to track the potential candidate genes associated with a particular trait of interest. The presence of high phenotypic variability in the nested synthetic introgression libraries coupled with the high marker density across the whole genome provided a strong base for the association mapping. Interestingly, the genes responsive to nutrient uptake under water stress (Diédhiou et al., 2008;Janicka-Russak and Kabała, 2015;Wang et al., 2017), shoot growth, root and plant development , nutrient uptake and transport of various nutrients (Takahashi et al., 2012;Wang et al., 2014b;Weng et al., 2020) reported to be collocated with 126 Mb genomic region on Chr 2D constituting 25 MTAs which stood out as hotspot for different traits (FRW, DRW, FSW, and DSW) in the present study. This indicates the positive interactions between root traits, nutrient uptake, and plant growth and development. The 7.7 Mb region on the short-arm of Chr 6A constituting 17 SNPs associated with GY showed collocation with the genes that were directly or indirectly involved in improving the GY in different cereal crops. These include the genes controlling flowering (Kania et al., 1997), panicle, and seed development (Jain et al., 2007;Li et al., 2011), GY (Terao et al., 2010), resistance to pathogenesis (Taniguchi et al., 2014;Wang et al., 2015;Niño et al., 2020) and abiotic stress tolerance (Brands and Ho, 2002;Palusa et al., 2007). The MTAs associated with different root traits, such as RSA, RV, tips, and forks in the present study were located near the earlier reported genes involved in regulating abscisic acid sensitivity and root growth development in Arabidopsis (Rodriguez et al., 2014) and adaptation under water stress conditions in wheat (Singh et al., 2017). Interestingly, the gene accelerating flowering in Arabidopsis (Hwang et al., 2019) was observed to be collocated with the SNP AX-95136655 associated with DTF on Chr 3B in the present study. The collocation of identified MTAs with earlier reported genes controlling the photosynthetic traits, root development, plant growth, nutrient uptake and transport, flowering, resistance to pathogenesis, and stress-responsive genes further confirms the contribution of these identified traits/MTAs in improving the N uptake/utilization and GY under N-limited conditions. The identified nitrogen irresponsive (NIR) breeding lines with favorable alleles in combination with the multiple traits might serve as potential donors for the development of N-efficient wheat varieties.

CONCLUSIONS
The nested synthetic introgression libraries covering extensive phenotypic variability coupled with huge genome coverage were used to identify the significant MTAs associated with nitrogen-use efficiency (NUE)-related traits in wheat. Significant phenotypic variations for the NUE-related traits, yield, and yield-related traits among genotypes, treatments, seasons, and their interactions (genotype x treatment, genotype x season, treatment x season, and genotype x treatment x season) were observed. Stable marker-trait associations (MTAs) identified for different traits measured in the present study comigrating with various genes associated with nitrogen (N) uptake/utilization and improving grain yield (GY) may help to harness their benefits in genomics-assisted breeding programs. The identification of N-efficient breeding lines may serve as novel donors in genomics-assisted introgression programs. The identification and introgression of superior haplotype improving NUE while maintaining GY using haplotype-based breeding may open new avenues in designing next-generation N-efficient high-yielding wheat varieties.

DATA AVAILABILITY STATEMENT
The original data and related information presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
NS and PC designed this study. AK provided the genotypic data of two populations and contributed to the development of nested introgression libraries. NS and MS conducted the field experiments. NS analyzed the data and wrote the manuscript. NS, SK, and PC provided resources. AK, MS, SK, VP-S, AS, AB, TB, and PC revised the manuscript. VP-S provided support with SPAD meter. All authors contributed to the article and approved the submitted version.