QTL analysis for nitrogen use efficiency in wheat (Triticum aestivum L.)

The genetic architecture of nitrogen use efficiency (NUE) and its two component traits i.e. NUpE (N uptake efficiency) and NUtE (N utilization efficiency) was studied using a bi-parental RIL mapping population derived from a cross HUW468 (high NUE)/C306 (low NUE). The mapping population, two parental genotypes and three check genotypes were evaluated under four different N levels (0, 60, 120, and 180 kg/ha) over three years. A genetic map containing 456 SNP markers (2571.38 cM length) was used for QTL analysis. Thirty six main effect QTLs (17 QTLs for NUE, 13 NUpE and 6 QTLs for NUtE) distributed on 12 chromosomes (1B, 1D, 2A, 2B, 3A, 4B, 5A, 5B, 5D, 6A, 6D, and 7A) were identified at 2.52–9.27 LOD scores. Individual QTLs explained 6.65–22.89% phenotypic variation. Multi-traits QTLs (Mt-QTLs) and epistatic QTLs involving first-order epistatic (QTL × QTL) interactions were also discovered. Candidate genes (CGs, as many as 737) were mined from QTL regions which were mainly involved in metabolic process, cellular process and catalytic activity, etc.; differential expression was observed for 49 CGs in roots and 34 in shoots. The CGs encoded important transcription factors, transporters, etc. having a role in NUE. QTLs and CGs reported in this study enriched the available knowledge. Seven QTLs (including three Mt-QTLs) and QTLs involved in six epistatic interactions are recommended for MAS for improvement of NUE in wheat.


Introduction
Wheat is one of the most important crop worldwide, and the wheat grain is being consumed daily by billions of people around the world. Although globally, wheat production has been able to meet the demand over several decades, low rainfall or droughts in major wheat producing countries have always been a major concern. Also, a decline in world-wide annual growth rate in production from 3% for several decades to the current > 1% has been a matter of concern (Ray et al. 2013). It is also estimated that in order to meet the demand of increasing population, globally, by the year 2050, wheat production should increase by at least 50%, despite the challenges of biotic and abiotic stresses (https:// www. opena ccess gover nment. org/ demand-for-wheat/ 83189/). In order to sustain and improve wheat production, besides other measures, one of the important measures is variability among wheat cultivars for response to nitrogen (N) fertilizer, both in terms of its uptake and utilization. Therefore, charactrization and utilization of the variation in N use efficiency (NUE) of future cultivars is an important step for maximizing production (yield) potential.
Like all other crops, wheat also needs a number of macro-(N, P, K, Ca, Mg, and S) and micronutrients (Cl, B, Fe, Mn, Cu, Zn, Ni, and Mb) for its optimum growth and development (White and Brown 2010). Among all these nutrients, N is the most critical for optimum plant growth and productivity (Razaq et al. 2017). The amount of N that is available in the soil in major wheat growing areas is generally a limiting factor for the optimum productivity of each individual crop including wheat (Elser et al. 2007). The extensive use of high doses of N fertilizers in almost all crops including wheat is also a matter of concern (see Balyan et al. 2016;Wang et al. 2021;Liu et al. 2022). Nitrogenous fertilizers influence directly or indirectly adequate food supply to half of the global population (Yadav et al. 2017).
Among cereals, three major staple crops including wheat, rice, and maize consume > 90.0% of the total nitrogenous fertilizers (Yadav et al. 2017). Estimates have shown that the global annual use of N fertilizers has been increasing over several decades and that it increased by 1.54% from 110.02 to 118.76 Mt during a period of recent five years (2015-2020) (http:// www. fao. org/3/ a-i6895e. pdf). However, the utilization of applied N is only 33% at the global scale with 42% in developed countries and 29% in developing countries (Raun and Johnson 1999;Ranjan et al. 2019; for reviews see Hawkesford 2017;Liu et al. 2022). A recent survey in China showed that in absolute terms, the N that remains unutilized (and lost) amounts to 52.5 kg N/ ha (range: 4.6-157.8 kg N/ha) (Tian et al. 2022). This unused soil N causes damage to the environment through surface water pollution, soil acidification, pollution of water bodies and increased emissions of greenhouse gas (GHG) (Yang et al. 2017;Tedone et al. 2018;Fradgley et al. 2021;Mălinaş et al. 2022). The application of high doses of N fertilizers also adds to the cost of production thus reducing the income of the farmers (Javed et al. 2022;Liu et al. 2022). Therefore, development of 'truly green' genotypes possessing higher NUE (grain production per unit of available N in the soil; Moll et al. 1982) will be desirable. In other words, this would mean higher yield under low/moderate-N inputs (Tian et al. 2016;Wang et al. 2021;Raigar et al. 2022;Liu et al. 2022). Improvement of NUE in crop cultivars is an important research area, not only for improvement of productivity of crop cultivars and reduction of the cost of production, but also for minimizing the use of fertilizer (without loss in productivity) to mitigate the adverse effects of fertilizers on the environment. Knowledge of the genetic architecture of NUE is a prerequisite for planning a strategy for developing these nitrogen use efficient cultivars.
It is known that NUE is a quantitative trait controlled by a large number of genes. However, relative to other traits, fewer studies on QTL analysis of NUE have been conducted in majority of crops including wheat, so that any study towards understanding the genetics of this trait should add to the knowledge. The present study is one such effort to improve our understanding about the genetics of this trait.
Genetic architecture of NUE and its two component traits has not been fully understood, although, some information is available. Using interval mapping (IM), 61 QTLs, and using genome-wide association studies (GWAS), 85 MTAs have already been reported (WheatQTLdb; Singh et al. 2021). However some of the studies conducted were based on data collected using hydroponics, making the utility of the results of a limited value for breeding (for details see Balyan et al. 2016;Ren et al. 2017Ren et al. , 2018Fan et al. 2018;Zhang et al. 2019;Brasier et al. 2020;Ranjan et al. 2021). Therefore, there is a need and scope for further studies for understanding the genetic architecture of NUE in Indian germplasm. Keeping this in view, the present study was planned, where a biparental population derived from a cross between two Indian wheat cultivars, namely HUW468 (high NUE) and C306 (low NUE) was used for QTL analysis. This material has never been utilized for a genetic study on NUE. Attempts have also been made for the identification of candidate genes (CGs) involved in pathways relating to NUE.

Plant material
A mapping population comprising 149 recombinant inbred lines (RILs) derived from a cross HUW468 (high NUE)/C306 (low NUE) was used for the QTL interval mapping. The seed material and SNP data of the mapping population were kindly provided by AK Joshi and VK Mishra of BHU, Varanasi, India. For developing RILs, cv. HUW468, a dwarf high yielding cultivar was used as the female parent and cv. C306, a tall and water stress tolerant genotype, was used as a male parent. The cv. HUW468 developed at BHU, Varanasi has the following pedigree: CPAN-1962/TONI// LIRA'S'/PRL'S' whereas, the cv. C306 was developed long ago by CCSHAU, Hisar and has the following pedigree: RGN/CSK3//2*C5 91/3/C217/ N14//C281.

Field experiments
The field experiments were conducted at the Indian Institute of Farming Systems Research (IIFSR), (29° 4′ N latitude and 77° 48′ E longitude; with an elevation of 237 m above MSL), at Modipuram, Meerut, India. A set of 159 genotypes [149 RILs, two parents (HUW468 and C306), five filler lines and three check genotypes (PBW343, HD2967, RAJ3765)] were evaluated using an augmented block design. Four different N levels (0, 60, 120, and 180 kgN/ ha) were used. The experiment was repeated for three consecutive Rabi seasons: 2016-2017 (E1), 2017-2018 (E2) and 2018-2019 (E3). Each experiment comprised 12 blocks, each block contained 16 genotypes (13 RILs/parental genotypes + three check genotypes). The 12th block in each experiment contained eight RILs, five fillers and three check genotypes. Each genotype was planted in 1.8 m 2 plots (four rows of 3 m with row to row distance of 0.2 m) with a seed rate of 125 kg/ha.

Fertilizer application
Nitrogen was applied in the form of urea (46% N) and included the following four treatments: (i) N0 = 0 KgN/ha (control), (ii) N60 = 60 kgN/ha, (ii) N120 = 120 kgN/ha, and (iv) N180 = 180 kgN/ha. Urea was applied in equal amounts at the following three developmental stages: (i) time of sowing, (ii) tillering, and (iii) anthesis. At the time of sowing, recommended doses of phosphate (in the form of single superphosphate) and potassium (muriate of potash or MOP) were also applied.

Recording of phenotypic data
The data in the field experiments were recorded and used for the estimation of the following three attributes following modifications of Moll et al. (1982) and Ortiz-Monasterio et al. (1997), made over by Haile et al. (2012).
where Ntf = Total N present in the plant (straw + grains) at maturity after treatment, Ntc = Total N present in the plant (straw + grains) at maturity of control treatment (0 kg ha −1 ), Ns = Nitrogen supplied.
where Gyf = grain yield of treatment with N fertilizer, Gyc = grain yield of control, and Ns = nitrogen supplied.
Straw N yield in kg/ha (SNY) was calculated as follows: SNY = Straw N% × raw yield (kg∕ha) where straw N% was estimated by Kjeldahl method (AOAC 1970) using the following steps: (i) Straw samples were ground and then digested by boiling in concentrated sulphuric acid. A transparent solution was obtained which contained ammonium sulphate. (ii) Ammonium was converted to ammonia by adding excess base to the ammonium sulphate solution.
(iii) Ammonia absorbed in the boric acid containing mixed indicator was determined by titrating it with a standard acid. The endpoint was indicated by the appearance of pink color.
where GPC = grain protein content, and GY = grain yield.
GPC (%) was estimated using grain analyzer-FOSS Grain Infratec 1241 (Graintec Scientific, Denmark). For this purpose, 100 g seed of each genotype/ RIL was used and mean GPC (%) for three replications per genotype/RIL was recorded. The measured GPC (%) was adjusted to 12% moisture content using the following formula: where M.C. = moisture content.

Genotypic data
Genotyping data for SNPs based on genotyping-bysequencing (GBS) for the RIL mapping population was kindly provided by Dr. Uttam Kumar, Borlaug Institute for South Asia (BISA), Ladhowal, Punjab, India.

ANOVA and descriptive statistics
ANOVA was conducted using R Studio. The fixed effect model was defined using the lm function of the stats package in R, where type-II ANOVA was generated using the Anova function available in the car package in R. Analysis was performed in the individual N levels over each of the three years and the pooled data (BLUP values). Descriptive statistics including frequency distribution, mean, range, standard error, coefficient of variation (CV%), and Pearson's correlation coefficients involving NUE, NUpE, and NUtE were obtained using SPSS 17.0 software. Best linear unbiased prediction (BLUP) value (i.e. pooled values of the three years data) for each of the three traits was calculated using package lme4 in R (https:// cran.r-proje ct. org/ web/ packa ges/ lme4/ index. html).

Construction of linkage map
Framework linkage map of the RIL mapping population was constructed using MultiPoint v.4.2. The genetic distances were estimated using Kosambi mapping function. MultiPoint first calculates the recombination frequency (RF) between each marker pair using maximum likelihood. Linked markers are then grouped together using the "bound together" function, with one marker selected as a representative marker of the group. This is followed by Vol.: (0123456789) clustering of markers into linkage groups by establishing a limit of RF followed by ordering the markers of each group.

QTL analysis
Single-locus QTL analysis was conducted by Composite Interval Mapping (CIM) approach using QTL Cartographer's Zmap QTL, Model 6 with a window size 10 cM (Basten et al. 2002). The Kosambi mapping function was used to change the recombination fraction into the genetic distance. The relative contribution of genetic component (R 2 ) was calculated as the proportion of the phenotypic variation explained (PVE). QTL analysis was conducted both in the individual N levels over each of the three years and the pooled data (BLUP values) of individual N levels over the years. In addition, multiple-trait analysis involving multi-trait composite interval mapping (Mt-CIM) and joint Mt-CIM was carried out for the detection of multi-trait QTLs (Mt-QTLs) using the module Zmapqtl available in QTL Cartographer. The QTLs with PVE% > 15% were considered as major QTL.
Two-locus analysis for first order epistatic interactions was conducted using software ICIM V 4.1.0.0 (Meng et al. 2015). For detection of epistatic QTL (E-QTLs), the probability in stepwise regression was set at 0.0001 and the scanning step was 5 cM. For each QTL, the position with the highest LOD score was considered as the most likely position of the QTL.
A threshold LOD score of 2.5 was used for declaring significant QTLs. The LOD score of 2.5 was chosen as a compromise threshold LOD to avoid too many false positives at a LOD score of 2.0 used in several earlier studies and too many false negatives at the thresholds calculated on the basis of permutations (Churchill and Doerge 1994).

Comparison of QTLs with previously reported QTLs/ MTAs/MQTLs
The most promising wheat QTLs were compared with known QTLs, MTAs, and meta-QTLs (MQTLs) reported in earlier studies. Each known QTL/MTA/ MQTL overlapping with a QTL identified in the present study was accepted to be co-located with QTL.

Identification of candidate genes (CGs)
Identification of CGs involved the following steps: (i) nucleotide sequences of markers flanking a QTL were utilized to identify the physical positions of QTL. For this purpose, blast (E-value ≥ 1E−100 and identity ≥ 95%) analysis of the nucleotide sequence was carried out against wheat reference genomic sequence (Triticum_aestivum_IWGSC_ensembl_release49) available on EnsemblPlants database. If the sequence of a marker was absent, the sequence of the adjacent genetic marker was used. (ii) Physical interval of each QTL region was calculated based on the genetic interval of the QTL region. The physical interval (in Mb, calculated from the coordinate information of the QTL) was divided by the genetic interval (in cM) and the distance in units of bases per cM was calculated. (iii) 1 Mb region on either side of the QTL (total 2 Mb interval) was used for identification of the candidate genes (CGs) underlying individual QTL regions.

Gene ontology (GO) and in silico expression analysis of CGs
GO analysis was conducted using Biomart tool available in Ensemble Plants. In silico expression of CGs was conducted using expression data available in Genevestigator software. The available expression data used in the present study belonged to Sharma et al. (2020). In this study four different wheat genotypes, namely OK1059060, BigSky, OCW00S063S, and NF97117 were evaluated at four N doses in the greenhouse.
The data on expression in Genevestigator was available as log2 transformed tpm (transcripts per million) value. Only genes (CGs) showing fold change (FC) ≥ 2 (up-regulation) or FC ≤ − 2 (downregulation) were accepted as showing significant differential expression in the form of fold changes estimated by comparing tpm values under low vs. sufficient N. The results of such differentially expressed CGs were depicted in the form of heatmaps generated using TBtools v1.0098726 software (https:// bio. tools/ tbtoo ls).

Phenotypic variation
Descriptive statistics including mean, range, SD and coefficient of variation (CV%) under individual N levels for each of the three years is given in the Supplementary Table S1; statistics on BLUP data are summarized in Table 1. In brief, the parental cv.
HUW468 was superior over the other parent C306 for all the three traits (NUE, NUpE, and NUtE), each following a decreasing trend in response to increasing N levels in both the parental genotypes. The RIL population also followed a similar trend (Table 1, Supplementary Table S1; Fig. 1). A wide range for each of the three traits was observed under all the three levels of N; CV(%) ranged from 10 to 20%. Violin plots of RIL population for NUE and its two component traits in individual years are given in Supplementary Fig. S1. The data for each of the three traits based on the BLUP values over three years showed continuous and normal distribution for all the three traits ( Fig. 1). Transgressive segregation was observed in both the directions for all the three traits under different levels of N.

Analysis of variance (ANOVA)
Results of ANOVA based on data for three years and three N levels are presented in Table 2 (Tables  for ANOVA involving each of the three individual  years are presented in Supplementary Table S2). As can be seen in Table 2, variations due to years, N levels and the genotypes for all the three traits were significant. The genotype × year interactions were also significant for all the three traits, but variation due to year × N level interaction was significant for only NUtE (Table 2). In the individual years, significant differences among the genotypes were observed for all the three traits under the three N levels (Supplementary Table S2).

Correlation among NUE, NUpE and NUtE
Values of Pearson's correlation coefficients among NUE, NUpE, and NUtE based on BLUP data of three years under three N levels are presented in Table 3. The three traits were positively correlated under three N levels. The correlation between the NUE and NUtE though positive but were only rarely significant. However, the correlations involving the NUE and NUpE were positive, significant and ranged from moderate to high.

Genetic linkage map
Out of 5717 polymorphic SNPs, only 456 SNP marker loci could be mapped on all the 21 linkage groups. The total map length was 2571.38 cM (Supplementary Table S3; Fig. 2). The sub-genome A was  Fig. 3). These included the following: (i) Seventeen (17) main effect QTLs for NUE, distributed on 10 chromosomes (1B, 1D, 2B, 3A, 5A, 5B, 5D, 6A, 6D, and 7A). Five of these QTLs were such which were detected at more than one N level in individual years and also BLUP data (   ccsu-3A.1 QTL was identified at N180; the highest PVE due to these QTLs was 22.9 and 17.1%, respectively and thus these may be considered as major QTLs. The remaining QTLs were moderate with PVE ranging from 6.7 to 13.39%. (iii) Six (6) main effect QTLs for NUtE, on chromosomes 2A, 2B, 5A, and 7A with LOD scores ranging from 2.5 to 3.8, included three QTLs on chromosome 5A alone (Table 4). Two of these QTLs were identified at more than one N levels including BLUP data. The remaining four QTLs were discovered at a single N level. All the QTLs were moderate with PVE ranging from 7.5 to 11.7%.

Co-located main effect QTLs for NUE and its two components
QTLs mapped on the same position of an individual chromosome are described by us as co-located, although some of them may also be detected using multi-trait QTL (Mt-QTL) analysis (see later). Over all, seven QTLs for NUE were co-located with QTLs for NUpE; these were located on six different chromosomes (1B, 3A, 5A, 5D, 6A and 6D); similarly, two QTLs for NUE were co-located with QTLs for NUtE, which were located on chromosome 5A.

CGs and GO terms
The genomic regions associated with main effect QTLs carried 737 CGs. The most important GO terms included the following: (i) biological processes: metabolic process, cellular process, response to stimulus, and biological regulation; (ii) molecular functions: catalytic activity, DNA binding and molecular function regulation etc. (Fig. 5   Vol.: (0123456789)

Discussion
The biparental population used in the present study exhibited sufficient variability for NUE, NUpE and NUtE, as shown by ANOVA (Table 2). Also, the phenotypic data exhibited normal distribution with transgressive segregation, suggesting that the population was suitable for the study, as also shown by the statistical analysis of phenotypic data collected on NUE and its two component traits (Fig. 1). Transgressive segregation also suggest that desirable QTLs were available in both the parents of mapping population. Since the three traits used in the present study exhibited positive correlations among themselves, multi-trait composite interval mapping (Mt-CIM) and joint Mt-CIM were also conducted to select Mt-QTLs, which can be used for simultaneous improvement of more than one correlated traits. In earlier studies also, simultaneous improvement of all the three traits has been recommended, although it will partly depend on the available variation for the three traits (Barraclough et al. 2010;Le Gouis et al. 2000;Brasier et al. 2020). The results of the present study were compared with the results reported in several earlier studies Habash et al. 2007;Guo et al. 2012;Xu et al. 2014;Zhang et al. 2019;Brasier et al. 2020;Ranjan et al. 2021) and the available information will be utilized for the discussion of the results of the present study (see later) related to QTLs for NUE and its component traits.

Main effect QTLs
For identification of QTLs, a LOD score of 2.5 was used as a compromise threshold, as explained in matreial and methods. Initially we did use a threshold LOD based on the concept of permutation (Churchill and Doerge 1994;Doerge and Churchill 1996), but in the present study, threshold LOD based on permuation differed for different traits and almost no QTLs were available, when threshold LOD based on permutation was used (not even in the genomic regions, where QTLs were reported in earlier studies on interval mapping; results not included). This was considered to be unrealistic and was presumably due to false negatives. The issue of calculation of threshold LOD has been widely discussed, and methods have been suggested for calculating threshold LOD scores (Lander and Botstein 1989;Van Ooijen 1999;Chen and Storey 2006). This issue of false positives and false negatives has also been widely discussed for QTL analysis involving genome-wide association studies (GWAS), where the recommended Bonferroni correction also gives a large number of false negatives (Perneger 1998;Kaler and Purcell 2019). Under these circumstances, a threshold LOD of 2.5 was used as a compromise to avoid both false positives and false negatives.
In the presenst study, QTLs were identified on 12 chromosomes. In earlier studies, QTLs have been reported on all the 21 chromosomes Laperche et al. 2006;Guo et al. 2012;Xu et al. 2014;Monostori et al. 2017;Hitz et al. 2017;Zhang et al. 2019;Ranjan et al. 2021). QTLs for other N-responsive agronomic and physiological traits have also been studied widely in the previous studies (for details see Balyan et al. 2016;Ren et al. 2017Ren et al. , 2018Fan et al. 2018;Ranjan et al. 2021); QTLs for these traits have also been reported on all the 21 wheat chromosomes (see Brasier et al. 2020). In earlier studies, NUE has also been estimated using root traits, shoot traits, dry matter, grain yield, reproductive/agronomic traits, N% in dry matter, etc. (wheatQTLdb; Singh et al. 2021Singh et al. , 2022. QTLs for dry matter have also been considered relevant for NUpE and those for photosynthesis may be relevant for NUtE. A comparison of QTLs detected in each of the three years also shows that none of the QTL was identified in more than one year. However, only one third of the QTLs (12/36) were identified in more than one N levels, suggesting the presence of QTL × N interaction in twothird (24/36) of the QTLs (Table 4). This is in agreement with earlier results, where no QTLs were detected in more than one environment Zhang et al. 2019;Xu et al. 2014;Brasier et al. 2020). These Q × E (including Q × N) interactions are important, making the utilization of these QTLs difficult. This also suggests that additional efforts are needed to identify stable QTLs with minimum QTL × E interactions.
A comparison of QTLs for three traits (17 for NUE, 13 for NUpE and 6 for NUtE) among themselves in our study confirms the earlier reported fact that maximum number of QTLs generally occur for NUE followed by QTLs for NUpE and NUtE. From earlier studies, a total of 87 QTLs were available for NUE, followed by 79 QTLs for NUtE and 33 QTLs for NUpE Laperche et al. 2006;Guo et al. 2012;Zhang et al. 2019;Brasier et al. 2020;Ranjan et al. 2021).
Some of the main effect QTLs detected in the present study for the three traits were located in the same genomic regions, where QTLs were also reported in earlier studies. For instance, four QTLs for NUE were co-located within the same genomic regions that were defined by QTL analysis in the two previous studies Brasier et al. 2020); these four QTLs (three on 5A and fourth on 7A) included the following: Qnue.ccsu-5A.1, Qnue. ccsu-5A.2, Qnue.ccsu-5A.3, Qnue.ccsu-7A.1). Of these, QTL Qnue.ccsu-5A.3 was also co-located with meta-QTL6 reported by Quraishi et al (2011), although none of the QTLs reported by us were colocated with any of the four meta-QTLs for NUE and its components recently reported by Saini et al. (2021). Thus, a majority of the main effect QTLs including the remaining 13 QTLs for NUE and all the QTLs for NUpE and NUtE reported by us were novel. A comparison of the QTLs reported in the present study with the MTAs reported in four earlier GWA studies for NUE and its components in wheat (Cormier et al. 2014;Hitz et al. 2017;Monostori et al. 2017;Xiong et al. 2019) revealed that two QTLs for NUE (Qnue.ccsu-1B.3 and Qnue.ccsu-7A.1) and a solitary QTL for NUpE (Qnupe. ccsu-3A.3) were also reported among the MTAs reported in the above GWA studies.
Among the three levels of N used in the present study, maximum number of QTLs for NUE and its two component traits were found at N120 followed by N180 and N60. However, with a view to improve NUE at the low/moderate soil N levels, the QTLs detected at more than one level of nitrogen having high PVE% are important for MAS; such QTLs need attention and include the following: (i) Two major main effect QTLs including Qnue.ccsu-1B.2 (PVE% = 6.7-19.6%) and Qnue.ccsu-5A.2 (PVE = 10.2-15.4%) for NUE discovered at more than one levels of N (N60/N120/N180) are important for improving NUE at different levels of available nitrogen. (ii) Another two QTLs for NUE, namely Qnue.ccsu-5D (detected in all the three N levels) and Qnue.ccsu-5A.1 (detected at N60 and N120), each having ~ 10% PVE are also important. (iii) A QTL for NUE (QTLs Qnue.ccsu-6D) having a PVE of 12.9% detected at N120 was co-located with a QTL for NUpE (Qnupe.ccsu-6D) that was detected at BLUP N60 and BLUP N120 is also important. (iv) Two QTLs for NUtE (Qnute.ccsu-5A.2 and Qnute. ccsu-5A.3) both located on chromosome 5A are also important (Table 4).

Co-located QTLs and Mt-QTLs
Each of the nine co-located QTLs (7 for NUE + NUpE; 2 for NUE + NUtE) identified during the present study may represent either a pair of two closely linked QTLs or a single QTL affecting more than one traits. Six such co-located QTLs include the QTLs, which were also listed as important in the above section on main effect QTLs, thus reaffirming their importance for breeding using MAS.
Seven Mt-QTLs that are listed in Table 5 are also important, since each such QTL controls more than one trait. Although, often in the published literature, such Mt-QTLs have been described as pleiotropic QTL (Fan et al. 2019), we prefer not to use this terminology, because of the following two reasons: First, because a QTL generally contains several genes, and second because, a locus unless Page 17 of 22 9 Vol.: (0123456789) ested for pleiotropy, may represent more than one closely linked genes. Interestingly, one of the Mt-QTL (WSNP3401 49927163 -WSNP1974 59193679 ) on 1B controlled all the three traits; the marker intervals in case of two Mt-QTLs on chromosomes 5A and 5D, one each was co-located with QTLs for the NUE and NUpE (Table 4). A novel and important Mt-QTL for NUE and NUpE (WSNP2007 46228336 -WSNP2809 46660512 ) on chromosomes 1B was detected at all the three N levels ( Table 5), suggesting that this may be used irrespective of the N-level used. Interestingly, the desirable alleles for all the above six Mt-QTLs were contributed by the superior parent HUW-468.

Epistatic interactions
In addition to the above main effect QTLs, 51 epistatic QTLs (E-QTLs) involved in 28 first order additive × additive interactions were also identified, which included 6 E-QTLs for NUE, 10 for NUpE and 12 for NUtE. These epistatic interactions may prove useful for wheat breeding, as also shown in several earlier studies involving many traits including NUE (Li et al. 2014;Ranjan et al. 2021). The fact that all E-QTLs (barring one exception on 2A, Qnue.ccsu-2A.4) did not include any main effect QTLs and the PVE due to individual interaction could be as high as 20.52%, suggests that epistatis interaction makes a significant contribution to the overall variability for the three traits under study (Supplementary Table S4). Also the fact that, out of 28 interactions,10 interactions belong to NUpE and 12 interactions belong to NUtE suggest that epistatic interactions are more important for NUpE and NUtE than NUE, for which only six interactions were available. However, the epistatic interactions, which have a PVE > 15% and are detected at N60 and N120 are more important for improvement of NUE under low/moderate N level and thus deserve further discussion. The following three interactions were most important for NUE: ccsu-7A.3 (N60 and N120) and Qnue.ccsu-5A.6 × Qnue.ccsu-5A.7 (BLUP N60). The first two of these three interactions were also responsible for NUpE at N120. Two other interactions were such which were only detected for NUpE at N60 and BLUP N120 and included the following: Qnupe.ccsu-5A.3 × Qnupe.ccsu-5B.1 and Qnupe.ccsu-5A.4 × Qnupe.ccsu-5B.2. Similarly, an important epistatic interaction Qnute.ccsu-1B × Qnute.ccsu-5A.4 was also detected for NUtE at all the three N levels. Epistatic interactions for NUE (assessed through shoot and root traits) were recently reported in wheat by Ranjan et al. (2021). Although, apparently, no example of the exploitation of epistasis for crop improvement through marker-assisted selection (MAS) is reported. However, importance of epistasis in MAS has been recommended and demonstrated through simulation by the group led by Jun Zhu in China (Liu et al. 2003). In this simulation study, it was demonstrated that MAS generally yields longer persistence response than that based exclusively on additive effect alone. It was also shown that neglecting epistasis could result in considerable loss in response, which becomes more pronounced at later generations. It has also been shown that incorporation of knowledge about epistasis can improve predictive ability for genomic selection in wheat (Raffo et al. 2022). We feel that in future additive × additive epistatic interactions will certainly be utilized to improve the outcome of a wheat breeding programme.

CGs and their in silico expression
Among a large number CGs identified during the present study, those encoding transcription factors (TF) and transporters are the most important, because particularly each of these TF in turn influence expression of multitude of genes including those involved in N uptake, transport and metabolism (see Raigar et al. 2022;Liu et al. 2022). The TFs encoded by CGs include a large number of domains listed in this article (Supplementary Table S5).
Our in silico analysis in the present study indicated that more genes (8-17) showed up-regulation in roots than the genes (3-4) showing upregulation in shoots of the four wheat genotypes, namely OK1059060, BigSky, OCW00S063S, and NF97117 (with higher NUE) under LN vs. SN inputs (Sharma et al. 2020). Some of these DE genes that encode TFs have relevance to NUE. (i) TF LBD (lateral organ boundary) that has contrasting expression in root and shoots of two genotypes (OK1059060 and NF97117) may mean that its role differ in relation to NUE in the two genotypes. In the past, however, the LBD (LBD37, LBD38 and LBD39) TF has been shown as negative regulators of NO 3 uptake and assimilation in Arabidopsis (Rubin et al. 2009;Zhu et al. 2022; for review see Liu et al. 2022) indicating its 9 Page 18 of 22 Vol:. (1234567890) role in NUE. (ii) The NAC TF, that plays a role in N remobilization in wheat, showed upregulation in shoot of OK1059060 under LN level (Zhao et al. 2015). The overexpression of TaNAC2-5A positively regulate the expression of multiple NPF (nitrate transporter 1/peptide transporter family) or NRT2 (nitrate transporter 2) and GS (glutamine synthetase) genes and increases NO 3 uptake, grain N concentration, N harvest index etc. (He et al. 2015). Thus the gene encoding the NAC TF is a useful resource for breeding crops with higher NUE.
Some other TF which did not show DE during our in silico study are also known to play a role in NUE. (i) The TFs belonging to different families (AUX/IAA, AP2/ERF, and WRKY) that are regulated by perception of NO 3 − showed higher expression at high/low N than at nil N in wheat (Effah et al. 2022;see Liu et al. 2022). (ii) Several genes that encode following TFs: GATA, ACT and bHLH. GATA TFs also have role in modulating N assimilation and metabolism balance between N and C (Hudson et al. 2011). The ACT TF can phosphorylate and inactivate AMT1;2, an ammonium transporter, under ammonium-sufficient conditions (Hsieh et al. 2018). Similarly, the over-expression of TabHLH1 exhibited drastically improved tolerance to N starvation in wheat (Yang et al. 2017). (iii) The TF NLP which contain PB1 domain regulates tissue-specific expression of genes involved in NUE and their expression is increased by NO 3 supply as well as N deficiency (Yu et al 2016;Kumar et al. 2018a, b;Jagadhesan et al. 2020). In our earlier study in wheat, it was shown that NLP genes are significantly up-regulated in the roots and shoots during N-stress (Kumar et al. 2018a, b). Since the expression of many of these TFs is also regulated by micro-RNA under N starvation (for a review see Zuluaga and Sonnente 2019), it is worth examining the role of the microRNA in NUE.
Among the different genes encoding transporters that were identified during the present study the SWEET sugar transporter showed differential expression in shoots/roots under low N level suggesting the possibility of the role of SWEET sugar transporters in NUE in wheat. In fact, in an earlier study in Arabidopsis, SWEET16 was shown to improve NUE (Klemens et al. 2013). Some genes encoding other transporters such as ABC transporter-like, Membrane transport protein, SLC26A/SulP transporter, auxin efflux carrier, plant type, and small auxin-up RNA did not show DE in our in silico expression analysis. Previously, ABC transporters showed upregulation in different tissues in the N-stressed roots in durum wheat suggesting their role in plant development and nutrition, organ growth etc. (Kang et al. 2011;Cai et al. 2012). However, DE CGs for nitrate and ammonium transporters so important for NUE were not discovered during the present study. The functional role of the remaining genes that were found to show DE during our in silico expression analysis need to be examined in future studies in relation to NUE in wheat.

Conclusions
The present study in wheat is one of the few studies involving identification of QTLs for NUE and its two component traits (NUpE and NUtE). Maximum number of main effect QTLs were detected for the NUE followed by QTLs for NUpE and NUtE. Several QTLs for NUE were co-located with QTLs for NUpE and rarely with QTLs for NUtE. Some of these co-located main effect QTLs were also multi-trait QTLs (Mt-QTLs). Nearly two dozen digenic epistatic interactions were also detected suggesting a major role of epistasis in the genetic control of NUpE and NUtE. The exploitation of the QTLs mainly detected at N60/N120 was recommended for use in MARS for improvement of NUE in wheat. Six important first order epistatic interactions were shortlisted for MAS for fixation of additive × additive interaction effects in early segregating generations. Important CGs associated with the detected QTLs that encode for TF and transporters play important role in NUE were also reported. There is also a need for the functional characterization of some of the differentially expressed genes detected during the present study. The information generated in this study improves our understanding of the genetic control of the complex traits such as NUE and its two components and is useful for the improvement of NUE through MAS in wheat. It is known that the higher grain yield and NUE in GA-insensitive dwarf rice are achieved if the ratio between the rice growth regulating factor-4 (GRF4) TF and the growth inhibitor DELLA protein is changed in favour of GRF4 abundance (Li et al. 2018). In future, the above ratio may be examined in the GAinsensitive dwarf wheats like HUW468 (higher NUE) or other similar genotypes to examine the possibility of a novel breeding strategies for further improvement of NUE in dwarf wheat. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.