Identification of Quantitative Trait Loci Associated with Nutrient Use Efficiency Traits, Using SNP Markers in an Early Backcross Population of Rice (Oryza sativa L.)

The development of rice cultivars with nutrient use efficiency (NuUE) is highly crucial for sustaining global rice production in Asia and Africa. However, this requires a better understanding of the genetics of NuUE-related traits and their relationship to grain yield. In this study, simultaneous efforts were made to develop nutrient use efficient rice cultivars and to map quantitative trait loci (QTLs) governing NuUE-related traits in rice. A total of 230 BC1F5 introgression lines (ILs) were developed from a single early backcross population involving Weed Tolerant Rice 1, as the recipient parent, and Hao-an-nong, as the donor parent. The ILs were cultivated in field conditions with a different combination of fertilizer schedule under six nutrient conditions: minus nitrogen (–N), minus phosphorus (–P), (–NP), minus nitrogen phosphorus and potassium (–NPK), 75% of recommended nitrogen (75N), and NPK. Analysis of variance revealed that significant differences (p < 0.01) were noted among ILs and treatments for all traits. A high-density linkage map was constructed by using 704 high-quality single nucleotide polymorphism (SNP) markers. A total of 49 main-effect QTLs were identified on all chromosomes, except on chromosome 7, 11 and 12, which are showing 20.25% to 34.68% of phenotypic variation. With further analysis of these QTLs, we refined them to four top hotspot QTLs (QTL harbor-I to IV) located on chromosomes 3, 5, 9, and 11. However, we identified four novel putative QTLs for agronomic efficiency (AE) and 22 QTLs for partial factor productivity (PFP) under –P and 75N conditions. These interval regions of QTLs, several transporters and genes are located that were involved in nutrient uptake from soil to plant organs and tolerance to biotic and abiotic stresses. Further, the validation of these potential QTLs, genes may provide remarkable value for marker-aided selection and pyramiding of multiple QTLs, which would provide supporting evidence for the enhancement of grain yield and cloning of NuUE tolerance-responsive genes in rice.


Introduction
Rice (Oryza sativa L.) is one of the most prominent staple food crops in the world and it has significantly contributed to global food security [1]. By 2050, rice production has to be increased by more than 60% to meet the rapid increase in food demand [2,3]. Approximately 90% of rice cultivation is carried out in irrigated and lowland rainfed ecosystems in Asia, where we still face problems such as decreasing arable land and freshwater availability, and increased labor costs and biotic and abiotic stress factors that impede production and productivity. Further, the drastic rise in fertilizer costs drives the search for suitable rice cultivars that are efficient in grain yield production [4][5][6][7][8]. Nutrient fertilizers are one of the most predominant factors that influence the genetic enhancement of grain yield productivity in most agricultural regions [6].
In the last four decades, the amount of nitrogen (N) applied to crops has risen 12 to 104 teragrams per year (Tg year −1 ), and 30%-50% of the N has been harvested in the grain and the remaining 70%-50% of the fertilizer has been lost through a combination of leaching, denitrification, volatilization, surface runoff, and microbial consumption [9]. Therefore, the excessive amount and long-term use of fertilizers cause severe environmental pollution, and this is one of the significant costly inputs for poor farmers [10][11][12][13][14]. Hence, the identification of crops with higher NuUE that are less dependent on NPK fertilizers (nitrogen, phosphorous, and potassium) is crucial for the sustainability of agriculture. Kant et al. [15] estimated that an increase of 1% in NuUE could save about US$ 1.1 billion annually. The Green Revolution period witnessed a modest increase of 13% in harvested rice area in Asia, whereas grain yield production more than doubled, from 240 to 513 million tons, with consumption of fertilizers high in NPK increasing from 6.7 to 69.0 million tons. According to the FAO electronic database, particularly in Asian countries, China, Sri Lanka, Indonesia, and India had the highest average fertilizer consumption of 565.2 kg ha −1 , 251.7 kg ha −1 , 185.141 kg ha −1 , and 157.48 kg ha −1 , respectively, from 2005 to 2014 (http://ricestat.irri.org:8080/wrsv3/entrypoint.htm#). Therefore, the identification of superior cultivars with improved nutrient use efficiency (NuUE) is an essential target research area for plant breeders, for attaining higher grain yield and also reducing production costs [16].
In the improvement of the tolerance of rice cultivars of single biotic and abiotic stresses, tremendous progress has been achieved [8,[17][18][19][20][21]. However, these tolerant varieties have not been able to succeed to attain higher grain yields under both the irrigated and rainfed ecosystems. To overcome these conditions, breeding aspects need to focus on developing multiple-stress-tolerant cultivars (MSTC) and improving nutrient use efficiency, which are vital for the sustainability of grain yield productivity, which could have a more significant impact on higher yield, especially under low-fertilizer-input conditions, besides being more beneficial to poor farmers. Several critical components of morphological, physiological, and agronomic traits are reducing grain yield and decreasing biomass content under nutrient-deficient conditions. These traits become altered due to changes in their molecular mechanism and physiological pathways, leading to their susceptibility [22]. The exploitation of rice genetic resources, using advanced genomic technologies, is essential for the identification of rice cultivars with NuUE for increasing crop grain yield under low-input conditions.
Developing rice varieties with stress tolerance and NuUE through conventional breeding approaches is extremely slow because of several factors that influence the molecular genetics and physiological mechanisms underlying low-input tolerance, the complexity of NuUE traits, and the absence of efficient breeding selection criteria [23][24][25]. The combination of advanced molecular marker-assisted breeding and conventional breeding platforms could speed up the breeding procedure for varietal development and identifying trait-associated QTLs. Over the past two decades, several types of traditional molecular markers, such as restriction fragment length polymorphism (RFLP), randomly amplified polymorphic DNA (RAPD), amplified fragment length polymorphism (AFLP) and simple sequence repeats (SSRs), have been used for QTL identification for NuUE-related traits under low-input conditions of nitrogen [26][27][28][29][30], phosphorus [7,[31][32][33][34][35][36][37][38], and potassium [39,40]. However, these molecular markers have some disadvantages vis-à-vis SNP markers, such as partial chromosome coverage, labor-intensiveness, low resolution, and higher cost. Among these various molecular markers, SNP markers have a comprehensive range of applications in the construction of genetic maps, uniform distribution, and cloning of QTLs. Their clear advantages, such as high density and accurate and reliable approaches, have been used to perform high-throughput genotyping based on SNPs across the genome [41][42][43][44]. In a genomics era, high-throughput genotyping with next-generation sequencing (NGS) and various array-based SNP detection platforms have become excellent tools for the dissecting of complex traits and identification of trait-associated genes and alleles in rice [45][46][47][48]. As compared to NGS, a SNP array can be used for many samples within a short period; it has a low cost and analysis of data interpretation is relatively easy [47,49]. The recently developed SNP array has been successfully used in diversity studies and genome-wide association studies (GWAS), and in the identification of numerous QTLs and genes in rice [43,[49][50][51].
Several QTLs may correspond to known genes in the N or P metabolic pathway, for example, Qyd-2b for N deficiency tolerance was located in the vicinity of the gene encoding cytosolic glutamine synthetase (GS1) [62], and Qyd-3b and Qpn-3 were nearby the genes for glutamate dehydrogenase (GDH2) [26]. Interestingly, Qyd-12 was detected only in low-P conditions and was co-localized with a major QTL, Pup1, on chromosome 12, which was reported to be involved in P absorption [63]. These results indicate that the QTLs specially detected under single N or P deficiency conditions may be involved in different pathways of N and P metabolism. The tightly linked markers have breeding potential in pyramiding elite genes and QTLs for N and P use efficiency.
Although N and P use efficient QTLs and linked traits were reported earlier in different genetic backgrounds of mapping populations, to date, there are no reports of different combinations of NPK fertilizers and their response to the foremost NuUE traits, such as grain yield response, agronomic efficiency, and partial factor productivity in rice. In the present study, taking advantage of the rapid development of SNP technologies, we used a 6K SNP array to genotype a BC 1 F 5 introgression line (IL) derived from a cross between Weed Tolerant Rice 1 (WTR-1), as the recipient parent, and Hao-an-nong (HAN), as the donor parent. In this study, we identified the NuUE-related QTLs by subjecting the population to varying nutrient rates and understanding the performance of yield and its components.

ANOVA and Interaction with the Environment
The fixed effect model was used for the analysis of variance with Satterthwaite denominator, and the summary of ANOVA is given in Table 2. Among the six nutrient conditions, significant genotypic effects were observed under -N (F = 1.44, p = 0.0005), -NP (F = 1.76, p < 0.0001), and -NPK (F = 2.47, p = 0.0007) conditions at the 1% level. In contrast, under 75N (F = 1.2, p = 0.053), -P (F = 1.38, p = 0.0658), and NPK (F = 1.07, p = 0.4293) conditions, the genotypes were not significantly different from each other. It was noted that the treatments showing significant genotypic effects are those that lack the nitrogen fertilizer component, indicating that some of these materials are nutrient use efficient. In the significant effect of fertilizers with genotype, the study used a −2 log-likelihood ratio test and found a significant environmental effect at the 1% level (χ 2 = 15.23, p = 0.0001). However, genotype and environment showed non-significant interactions (χ 2 = 0.00, p = 0.9992) ( Table 3). This suggests that the different fertilizer applications significantly affect yield performance of the ILs as a whole. On the other hand, the non-significant genotypes by environment interactions indicate the consistent performance of genotypes across environments.

Analysis of Agronomic Efficiency (AE) and Partial Factor Productivity (PFP)
In the NuUE study, agronomic use efficiency and partial factor productivity were calculated using the fertilizer application rate of NPK (160-50-50 kg ha −1 ). The computation using the equations (i to iii) revealed 15.06-34.11 kg of grain for each kg of N derived from an NPK fertilizer given in 3.2:1:1 ratio, respectively. However, using this equation, we found 117 ILs with >15 kg of grain per kg −1 N (NPK), along with 33 ILs (-P) and 161 ILs (75N). These were also compared with parents and we identified 42 ILs, 23 ILs, and 79 ILs that had a higher AE with respect to each equation (Table 4). Similarly, analysis of partial factor productivity helped to identify 16 ILs (NPK), 4 ILs (-P), and 151 ILs (75N) that had a >50 kg of grain per kg −1 N (NPK) condition. In comparison with parents, 25 ILs (NPK), 61 ILs (-P), and 117 ILs (75N) showed more than 50 kg of grain per kg −1 N in each condition.

Construction of Linkage Map and Segregation of SNP Markers
A total of 704 high-quality polymorphic SNP markers, with an average for each chromosome of 58 SNP markers, were used to genotype the 230 ILs, and a high-density genetic linkage map was developed covering 1526.8 cM, with an average of 127.2 cM per SNP marker ( Table 5). The logarithm of the odds (LOD) thresholds and explained phenotypic variance for six nutrient conditions ranged from 2.52% to 17.76% and 5.87% to 34.68%, respectively. Using single marker analysis (SMA), a total of 261 QTLs were mapped on all 12 chromosomes, with an average of 21 QTLs for yield attributed to seven key component traits of NuUE in rice. Among the 12 chromosomes, the highest number of QTLs was located on chromosome 6 (26 QTLs), whereas the lowest number of QTLs was observed on chromosome 12 (6 QTLs), and all QTLs are listed in Table S1. The hotspot QTL regions are shown in Figures 2 and 3. The majority of the QTLs were associated with 1000-Gwt, and PSPF traits had a negative additive effect of 61.9% and 61.5%, indicating that alleles from the recipient parent WTR-1 contributed to increasing phenotype.

Construction of Linkage Map and Segregation of SNP Markers
A total of 704 high-quality polymorphic SNP markers, with an average for each chromosome of 58 SNP markers, were used to genotype the 230 ILs, and a high-density genetic linkage map was developed covering 1526.8 cM, with an average of 127.2 cM per SNP marker ( Table 5). The logarithm of the odds (LOD) thresholds and explained phenotypic variance for six nutrient conditions ranged from 2.52% to 17.76% and 5.87% to 34.68%, respectively. Using single marker analysis (SMA), a total of 261 QTLs were mapped on all 12 chromosomes, with an average of 21 QTLs for yield attributed to seven key component traits of NuUE in rice. Among the 12 chromosomes, the highest number of QTLs was located on chromosome 6 (26 QTLs), whereas the lowest number of QTLs was observed on chromosome 12 (6 QTLs), and all QTLs are listed in Table S1. The hotspot QTL regions are shown in Figures 2 and 3. The majority of the QTLs were associated with 1000-Gwt, and PSPF traits had a negative additive effect of 61.9% and 61.5%, indicating that alleles from the recipient parent WTR-1 contributed to increasing phenotype. Chr-chromosome; Kb-kilo base pairs; cM-CentiMorgan. Out of the 261 QTLs identified, 62.5% (163 QTLs) had a positive additive effect contributed by WTR-1, and the remaining 37.5% followed with the negative allele from donor parent Hao-an-nong. QTLs with >20% phenotypic variation explained (PVE) were determined as major QTLs with an LOD score of >9.50, while others were considered as minor QTLs. In summary, a total of 261 QTLs for seven traits under six different combinations of NPK nutrients were identified, out of which 49 were major QTLs and 212 were minor QTLs (Table S1). Out of these NuUE conditions, a large number of QTLs were observed in -P (61 QTLs), -NP, and -NPK (50 QTLs), and followed with 75N Out of the 261 QTLs identified, 62.5% (163 QTLs) had a positive additive effect contributed by WTR-1, and the remaining 37.5% followed with the negative allele from donor parent Hao-an-nong. QTLs with >20% phenotypic variation explained (PVE) were determined as major QTLs with an LOD score of >9.50, while others were considered as minor QTLs. In summary, a total of 261 QTLs for seven traits under six different combinations of NPK nutrients were identified, out of which 49 were major QTLs and 212 were minor QTLs (Table S1). Out of these NuUE conditions, a large number of QTLs were observed in -P (61 QTLs), -NP, and -NPK (50 QTLs), and followed with 75N (42 QTLs) and -N (40 QTLs), and the fewest were detected in NPK conditions with 18 QTLs. The PVE of significant major QTLs ranged from 20.25% to 34.68% at 9.55 to 17.76 LOD value, whereas, for minor QTLs, the explained PV ranged from 5.87% to 19.98% at 2.52 to 9.49 LOD value. The distribution of all major and minor QTLs, represented on all 12 chromosomes, is shown in Figure 3, and the hotspot QTLs (>5 QTLs), identified over nine chromosomes (chr01, chr02, chr03, chr04, chr05, chr07, chr08, chr09, and chr11), were linked to six essential nutrient traits (PFP, GY, BY, FGP, 1000-Gwt, and PSPF) in rice.  Figure 3, and the hotspot QTLs (>5 QTLs), identified over nine chromosomes (chr01, chr02, chr03, chr04, chr05, chr07, chr08, chr09, and chr11), were linked to six essential nutrient traits (PFP, GY, BY, FGP, 1000-Gwt, and PSPF) in rice. Out of the 704 polymorphic markers, SNP_2_4481943 had the highest LOD value recorded (17.76), which explained PV of 34.68%, followed by two other markers, SNP_10_6149421 (LOD 16.28) and SNP_8_10073191 (LOD 15.31), which explained PV of 32.32% and 30.73%, respectively. Interestingly, these QTLs were significantly associated with 1000-Gwt under -NPK conditions, and their negative value of the allele was exhibited from donor parent Hao-an-nong. In contrast, the lowest LOD (2.52) was recorded in SNP_6_9836381, which explained PV of 5.87%, followed by three other SNP markers, SNP_12_14936674 (2.55), SNP_12_7445812 (2.55), and SNP_12_5851455 (2.53), which explained PV of 5.92%, 5.92%, and 5.89%, respectively. These four QTLs were linked with different key traits such as FGN, AE, BY, and GY under -NP, 75N, -N, and -NPK conditions. The lowest LOD value of two markers (SNP_6_9836381 and SNP_12_14936674) was contributed by the Out of the 704 polymorphic markers, SNP_2_4481943 had the highest LOD value recorded (17.76), which explained PV of 34.68%, followed by two other markers, SNP_10_6149421 (LOD 16.28) and SNP_8_10073191 (LOD 15.31), which explained PV of 32.32% and 30.73%, respectively. Interestingly, these QTLs were significantly associated with 1000-Gwt under -NPK conditions, and their negative value of the allele was exhibited from donor parent Hao-an-nong. In contrast, the lowest LOD (2.52) was recorded in SNP_6_9836381, which explained PV of 5.87%, followed by three other SNP markers, SNP_12_14936674 (2.55), SNP_12_7445812 (2.55), and SNP_12_5851455 (2.53), which explained PV of 5.92%, 5.92%, and 5.89%, respectively. These four QTLs were linked with different key traits such as FGN, AE, BY, and GY under -NP, 75N, -N, and -NPK conditions. The lowest LOD value of two markers (SNP_6_9836381 and SNP_12_14936674) was contributed by the negative allele from donor parent Hao-an-nong, and two other markers (SNP_12_7445812 and SNP_12_5851455) had a positive allele from recipient parent WTR-1.

Agronomic Efficiency (AE)
Agronomic efficiency (AE) was determined in the experiment on NuUE, and we found four QTLs (qAE_2.1, qAE_4.1, qAE_6.1 and qAE_12.1) that were located on chromosomes 2, 4, 6, and 12, respectively ( Table 6). For the agronomic efficiency of applied nitrogen in terms of -P conditions (applied nitrogen and potassium fertilizer), three QTLs were detected with an LOD score and phenotypic variation of 2.77, 4.01, and 4.52 and 6.43%, 9.17%, and 10.27%, respectively, with a positive additive effect of 3.16, 2.13, and 2.28, which indicates that the progeny carried the trait from the recipient parent WTR-1. QTL qAE_6.1 consisted of a peak marker (SNP_6_9977282) that was also associated with GY and PFP under -P conditions, with explained phenotypic variation of 17.6%. In 75N conditions, one QTL (qAE_12.1) was detected with a peak marker (SNP_12_14936674) located on 14936674 bp of chromosome 12, whereas LOD score and phenotypic variation were 2.55% and 5.92%, respectively, since the additive value of -2.8 shows that the trait in the progeny came from donor parent Hao-an-nong.

Partial Factor Productivity (PFP)
A total of 22 QTLs were identified across all chromosomes, except chromosome 12, under two different nutrient conditions (-P and 75N). Eleven QTLs were associated with PFP with LOD scores ranging from 3.39 (SNP_5_5588965) to 12.15 (SNP_10_6149421), which explained phenotypic variation from 7.81% to 25.28% in -P conditions. Among these 11 QTLs, the highest LOD value of three peak markers, SNP_2_4481943 (LOD 11.68), SNP_8_8437588 (LOD 9.9), and SNP_10_6149421 (LOD 12.15), was contributed by the Hao-an-nong allele on chromosomes 2, 8, and 10, and eight QTLs remained on chromosomes 1, 3, 4, 5, 6, 7, 9, and 11 that had a positive additive value, indicating the contribution of recipient parent WTR-1. Under -P conditions, one of the peak markers (SNP_1_20345712) was located at the 20345712 bp position on chromosome 1, which is associated with BY and GY, with an LOD score of 4.01 and 8.63, explaining a PV of 9.17% and 18.71%, respectively (Table 6). Similarly, we discovered 11 significant QTLs for PFP under 75N conditions. Of the 11 QTLs detected, four (qPFP_2.1, qPFP_5.2, qPFP_8.1, and qPFP_10.1) were contributed by the donor parent allele of Hao-an-nong, and the remaining seven QTLs by the WTR-1 allele. An average LOD value of 7.22 and PV of 15.82% were explained by 11 QTLs. The two highest peak markers were SNP_5_15469279 and SNP_2_4342883, located on chromosomes 2 and 5, respectively, which explained PV of 20.25% and 20.91%. In total, 52 QTLs were identified for six nutrient conditions and they were distributed on all chromosomes except for chromosome 12. These QTLs explained phenotypic variation ranging from 5.98% to 29.76%. Of these 52 QTLs, 6 QTLs for 75N, 3 QTLs for NPK, 10 QTLs for -P, and 11 QTLs for each nutrient condition (-N, -NP, and -NPK) were significantly expressed, with an LOD score ranging from 2.57 to 14.73. Under -N deficiency conditions, four QTLs (qPSPF_3.3, qPSPF_6.4, qPSPF_8.2, and qPSPF_9.1) and another three QTLs (qPSPF_3.4, qPSPF_6.5, and qPSPF_9.2) for -NPK conditions were recorded with PVE of >20%. The peak markers were located on chromosomes 3, 6, 8, and 9 with the favorable allele contribution from WTR-1.

Discussion
NuUE is a complex trait influenced by several factors, and it is considered a vital trait to improve rice grain yield productivity in marginal and rainfed lowland areas. Rice breeding programs need to incorporate NuUE traits as this helps resource-poor farmers to save on fertilizer, and allows maximization of valuable resources to increase profitability. To date, several morphological and agronomic traits have been identified, and they can be used as indicators of nitrogen (N) [52,[64][65][66], phosphorus (P) [7,38,[67][68][69][70], and potassium (K) [16,39,40,[71][72][73] deficiency tolerance under both field and hydroponic conditions. Using different genetic backgrounds of mapping populations, such as RILs, NILs, BILs, ILs, CSSLs, DHs, and BC 2 F 4 , there are several reports on QTL identification with independent studies of low-N, -P, and -K conditions. The previously identified QTLs were mapped on different chromosomes using low-density linkage maps that were constructed by PCR-based molecular markers, such as RFLP, AFLP, SSR, and STS [26,29,31,34,39,48,52,55,58,67,[74][75][76]. However, in this study, we used SNP markers to identify the QTLs.
For the identification of QTLs, single-marker analysis, simple interval mapping, and composite interval mapping are widely used methods [77,78]. The development of nutrient use efficient selective ILs, through an early generation backcross breeding approach, allows both the detection of QTLs for NuUE traits and the simultaneous development of promising breeding materials into varieties. This is the strength of selective introgression breeding first proposed by Tanksley [78], and later demonstrated for rice by Li et al. [79]. However, for QTL detection in selective introgression backcross breeding populations, we could use single marker analysis (SMA) and association mapping with single nucleotide polymorphism (SNP) markers, and use a higher threshold of LOD value of 2.5 to declare a strong QTL.

Analysis of Critical NuUE Traits
The phenotypic evaluation of 230 BC 1 F 5 ILs varied significantly in six different NPK nutrient combinations: N, -P, -NP, -NPK, 75N, and NPK. Grain yield is a quantitative trait and a highly complex character for all crops [80], and is the essence of any breeding program. Various morphological and physiological plant traits contribute to grain yield. Yield components are inter-related with each other, indicating a complex chain of relationships, which is highly influenced by the environment [81]. The breeding strategy in rice mainly depends on the degree of associated characters as well as the magnitude and nature of variation [82,83]. Based on a higher grain yield level, 18 ILs for NPK, 23 ILs for 75N, 4 ILs for -P, and 1 IL for both -NPK and -NP had more than 40 g/plant, and in -N conditions the lowest GY was observed (34.96 g/plant).
An understanding of the correlations between traits is of great importance in breeding programs, especially if the selection of one of them is impaired by low heritability or difficulties of measurement and identification [84]. Information on trait correlations has been helpful as a basis for selection in breeding programs. Ashura [85] showed a positive correlation between the number of filled grains per panicle, number of panicles per plant, and 1000-grain weight for grain yield. A correlation study enables breeders to understand the major traits for which selection can be based on population improvement. GY had a significant positive correlation with BY and FGN in NPK; FGN in 75N; GWT, BY, and FGN in -P; PSPF, BY, and FGN in -N; and BY and FGN in -NP (p < 0.01) conditions. A negative correlation was recorded for GY with 1000-Gwt in NPK and for BY in 75N conditions. Under all six nutrient conditions, GY showed a negative relationship with BY and the remaining traits, GWT, PSPF, and FGN, had a significant positive correlation between the traits. Highly correlated traits were indicated to have a common genetic basis, suggesting that these eight relative phenotypic traits could be used for the evaluation of P-deficiency tolerance in BILs at the seedling stage [86].

Promising Traits of AE-Associated QTLs
The application of AE and PFP traits was essential to plant nutrients for an optimum quantity and right proportion, through the correct method and time of application, and this is the key to increased and sustained crop production [87]. These two traits are important for the measurement of nutrient use efficiency, and this also provides an integrative index that quantifies total economic output relative to the use of all the nutrient resources in the system [88]. In 1996, Cassman et al. [10] defined PFP and AE and showed that they could be increased by raising the amount, uptake, and use of available nutrients, and by increasing the efficiency with which applied nutrients are taken up by the crop and used to produce grain. In this study, we identified 42 ILs with higher AE in NPK, 23 ILs in -P, and 79 ILs in 75N, as compared with their parents, and they were explicit in >15 kg grain kg −1 N conditions. Yoshida [89] and Cassman et al. [10], respectively, estimated AE to be 15-25 kg kg −1 and 15-20 kg kg −1 in the dry season in farmers' fields in the Philippines. In a similar way, Wen-Xia et al. [90] mentioned the AE in two kinds of rice: Jinzao had an AE range of 8.02-20.14 kg grain kg −1 N and Shanyou63 had an AE range of 3.40-18.37 kg grain kg −1 , differing with N management. In hybrid rice, AE for applied P was 5.2 kg grain kg −1 P and for applied K was 11.8 kg grain kg −1 K, where the fertilizer application rate for NPK was 200-75-200 and 200-150-200 kg ha −1 , respectively; AE for applied P in non-hybrid rice was 2.3 kg grain kg −1 P and 4.7 kg grain kg −1 P, where the fertilizer rate was the same [81].
Four novel QTLs were identified for AE (qAE_2.1, qAE_4.1, qAE_6.1, and qAE_12.1) on chromosomes 2, 4, 6, and 7, which explained phenotypic variation ranging from 5.93% to 10.27%. For these four AE QTLs, each SNP marker position was analyzed at both sides of the 500 kp (up-and down-stream of 1 Mb) regions on the chromosome (Table S3). With the respective positions of the 1 Mb region on chromosome 2, seven genes were identified influencing diverse functional roles as a defense to pathogen resistance and physiological mechanisms. Of the seven genes revealed, BiP3 for BLB resistance [91]; OsHPL3 for multiple biotic disease resistance, such as brown planthopper, striped stem borer, and BLB [92]; BiP1 (Os06g0622700) for increased seed storage protein, starch content in the endosperm, and also stress responses [93,94]; OsHPR1 (Os02g0104700) involved in photo-respiratory metabolism [95]; and OsNOA1 regulating chlorophyll biosynthesis, plastid development, and Rubisco formation in a temperature-dependent manner [95], were identified on chromosome 2. Similarly, on chromosome 6, genes Pi9 and Pi2 were noticed [96], which are related to blast disease resistance, along with another three QTLs (qPHw6-2, qSFWd6, amy6-1) related to drought tolerance and amylose content [97,98]. Interestingly, a major QTL for P deficiency tolerance that was previously mapped on chromosome 12 was found to be closely associated with peak marker SNP_12_14936674 [31], and in the same region of markers was linked with drought tolerance QTLs [99].
In PFP analysis, 16 ILs in NPK, 4 ILs in -P, and 151 ILs in 75N conditions yielded >50 kg grain kg −1 N and, in comparison with parents, 25 ILs (NPK), 61 ILs (-P), and 117 ILs (75N) yielded more than 50 kg grain kg −1 N in each condition. Amanullah et al. [100] showed, with a maize crop, a PFP of applied N of 36.62 kg grain kg −1 N and an AE of applied N of 22.49 kg grain kg −1 N. After applying DAP and SSP in fields, the AE of two fertilizer applications resulted in 13.01 and 13.71 kg grain kg −1 P, and PFP resulted in 63.58 and 61.92 kg grain kg −1 P. In this study, 22 QTLs for PFP traits were identified under 75N (11 QTLs) and -P (11 QTLs) conditions, with a PV ranging from 8.41% to 20.91% and 7.81% to 25.28%, respectively. Out of these 22 QTLs, six QTLs at SNP marker positions were analyzed on both sides of 500 kb (up-and down-stream of 1Mb) regions, which were explained by PV of more than 20% on chromosomes 2, 4, 5, 8, and 10. Five genes and six QTLs were identified in the 1 Mb region associated with PFP traits. On chromosome 2, three genes were identified as OsWRKY71 [101], Os4CL3 [102], and OsWRKY45 [103,104], which were functionally related to defense signaling molecules, such as SA, Me, JA and lignin biosynthesis. Two QTLs (qDSR_8 and qALRR_8) on chromosome 8 [105,106] and three QTLs (qRRE_10, qRFW_10, and qRCCL_10) [107,108] and one gene, OsAT1/Spl18 [109], located on chromosome 10 revealed having tolerance of aluminum and alkaline stress, and also showed blast disease and ultraviolet-B resistance in rice. On chromosome 5, gene OsWRKY45 was found to be closely associated with SNP_5_15469279, and it is involved in multiple biotic and abiotic stress tolerance/resistance, such as BLB, blast, sheath blight, drought, salinity, and cold [103,104,110]. A comprehensive literature survey revealed that a majority of the QTLs associated with low phosphorus were reported on chromosomes 1, 2, and 12 (recently reviewed by Mahender et al. [7], Ali et al. [20], and van de Wiel et al. [70], whereas, for low nitrogen conditions, the majority of morpho-physiological trait-linked QTLs in rice were located on chromosomes 3, 5, and 8 [34,[52][53][54]76]. Under six different nutrient conditions, 24.1% of the QTLs were associated with 1000-Gwt, followed by 19.9% for PSPF, 18.7% for GY, 15.3% for BY, 11.8% for FGP, 8.4% for PFP, and 1.5% for AE (Figure 3).

Consistency and Comparisons of Major QTLs across Different Genetic Backgrounds
In this study, a total of 261 QTLs were identified using a 6K SNP array-based genetic linkage map analysis, in 230 BC 1 F 5 introgression lines, that were tested under six different nutrient conditions. Out of 261 QTLs, 49 major QTLs were found with more than 20% PVE and LOD threshold value ranging from 9.44 to 17.76 distributed across all chromosomes, except for chromosomes 7, 11, and 12. Of these 49 QTLs, the highest number (27 QTLs) was associated with 1000-Gwt and the lowest number (6 QTLs) with PFP. For other traits, nine QTLs were associated with GY and seven QTLs with PSPF, which were detected across five nutrient conditions: 75N, -P, -N, -NP, and -NPK. Interestingly, out of these 261 QTLs, 14 QTLs were consistently detected under more than four nutrient treatments, indicating that these were critical traits that were controlled by similar genes under different nutrient conditions. The consistencies of QTLs were identified on nine chromosomes: 1, 2, 3, 5, 7, 8, 9, 10, and 11.
Interestingly, the identified associated QTLs governing seven traits (1000-Gwt, FGN, PSPF, BY, GY, AE, and PFP) on chromosome 6, under phosphorus-deficient conditions, significantly contributed to PVE ranging from 7.51% to 27.49%. Earlier independent studies of Ni et al. [31] and Wissuwa et al. [55] mapped major QTLs for RTA, RSDW, and RRDW on the same chromosome 6, and a group of P-responsive genes and transcription factors were also located in this region [113], which may confer tolerance of P deficiency [117]. However, in comparison to Pup1 on chromosome 12, major QTLs and P-responsive genes were mostly located on chromosome 6, which indicates that both have independent genes and regulatory pathways [118].
In the deficiency of -N, -P, -NP, and -NPK conditions, four QTLs for GY were identified in -P (qGY_2.4, qGY_4.2, qGY_8.4, and qGY_10.5), two QTLs in -NP (qGY_5.4 and qGY_10.4), and one QTL in -NPK (qGY_5.3), which were detected on different chromosomes (2, 4, 5, 8, and 10) from the analysis of 49 major QTLs. These QTLs were explained by their PV ranges from 20.25% to 25.28%. This was in contrast to Yue et al. [52], who reported three QTLs for 1000-Gwt on chromosome 3 and 7 and another three QTLs for GY located on chromosomes 1, 4, and 9, and the phenotypic variation explained by these QTLs ranged from 4.93% to 26.73% and 5.73% to 6.80%, respectively. However, in the present study, on chromosome 4, peak marker SNP_4_21833014 was found to be close to qGYP-4 and was also flanked by RM273-RM241 [52]. In the present study of the genomic region of SNP markers from SNP_1_195334 to SNP_1_23839187 on chromosome 1, five QTLs were identified (qBY, qGY, qFGN, qPSPF, and q1000Gwt), which explained PV ranging from 7.18% to 26.83% under N-deficient conditions. However, in the same genomic region of chromosome 1, several QTLs governing relative grain yield, grain weight, 1000-Gwt, and nitrate transporter, such as NRT 2.1 and OsNRT2.3,b are significantly involved in nitrogen deficiency tolerance and improving N uptake and enhancement of grain yield under nitrogen deficiency [119][120][121][122][123]

QTL Hotspots
Association with several traits within a single genomic region has immense potential value for the development of desired target trait enhancement through breeding and marker-assisted selection (MAS) applications. A single gene may affect more than two traits through pleiotropy or with closely linked genomic loci [124,125]. In the present study, 14 hotspot QTLs were identified to be located on chromosomes 1, 2, 3, 4, 5, 7, 8, 9, and 11, which had PVE ranges from 5.99% to 34.68% (Figure 3). Of these 14, we revealed the top four QTL harbor regions on chromosomes 3, 5, 9, and 11, and were designated as QTLs harbor-I to -IV. Each of these QTL harbor regions contained more than 10 QTLs on the same genomic regions of SNP markers, and significantly carried the positive allele from recipient parent WTR-1. The top four hotspot QTLs (a total of 45 QTLs) located on chromosomes 3, 5, 9, and 11 had PVE ranging from 6.09% to 26.97% and LOD values of 2.62 to 13.11. The four chromosomal (3, 5, 9, and 11) regions were shared by nitrogen (N) [34,52], phosphorus (P) [39,48,55,[59][60][61], and potassium (K) [39] deficiency QTLs reported in different mapping populations in rice.
Significant reports exist on chromosomes 3, 6, and 11 involved in phosphorus deficiency tolerance in rice. Using BILs derived from inter-specific crosses of Oryza sativa L. X O. rufipogon Griff. [38], F 3 lines from the crosses between P-deficiency-tolerant variety NERICA10 and sensitive variety Hitomebore [126], and 271 introgression lines [60], identified QTL clusters for SDW, BY, and P uptake on chromosomes 3 and 11, which showed high significance for P-deficiency tolerance. Comparison with previous studies showed a majority of the candidate genes for P-deficiency tolerance (Pup1) and PSTOL1 located on chromosome 12 [63,112,118,127] and P-responsive genes (such as OsPTF1) on chromosome 6 [117]. Our present study revealed that a major region of chromosomes 3, 5, 9, and 11 might provide novel loci to discover candidate genes related to P-deficiency tolerance, and these genes/alleles might play a co-regulated role for P-deficiency tolerance. These significant results would be more helpful to breeders and biotechnologists to understand the genetic and molecular basis of the physiological mechanisms of P-deficiency tolerance in rice.
According to the physical positions of one of the QTL hotspot regions, SNP marker SNP_3_853802 was associated with the highest number of 14 QTLs that were located on chromosome 3, and were distributed as five QTLs in -P (q1000gwt_3, qPFP_3, qBY_3, qFGN_3, and qGY_3), two QTLs in -N (q1000gwt_3 and qGY_3), two QTLs in -NP (qBY_3 and qGY_3), three QTLs in -NPK (q1000gwt_3, qBY_3, and qGY_3), and one QTL in 75N and NPK (q1000gwt_3). Hạnh et al. [141] identified three hotspot regions with response to low-nitrogen QTLs, which were flanked by RM265-RM165 on chromosome 1, RM3199-RM514 on chromosome 3, and RM080-RM281 on chromosome 8, which were significantly associated with the traits as total fresh weight of leaf blades (FW), RDW, SDW, nitrogen concentration in sheaths plus stem (NS), nitrogen concentration in leaf blades (NL), plant height (PH), chlorophyll content index (CCI), NUE, physiological nitrogen use efficiency (pNUE), and agNUE in RIL populations of IR64 and Azucena. In another study, Senthilvel et al. [58], using doubled haploid (DH) lines of IR64/Azucena in a pot experiment with three N doses (native, 0 kg/ha −1 ; normal, 100 kg/ha −1 ; and high, 200 kg/ha −1 ), identified seven main-effect QTLs that were associated with NUE traits and plant grain yield on chromosome 3 [58]. In similar studies, with two different N fertilizer rates (0 N and 130-135 kg N ha −1 ), in field conditions with RILs derived from two Oryza sativa-ssp. indica rice varieties (Zhenshan97/Minghui63), Wei et al. [76] reported a major QTL for grain yield (qRGY3) on chromosome 3, and another two QTLs (qRGY7 and qRGY11) on chromosomes 7 and 11, detected with PV of 10.8% in 2006 and 16.0% in 2007. Moreover, three QTLs (qRBM9-1, qRBM9-2, and qRBM10) for biomass yield on chromosomes 9 and 10 together explained 33.6% of the total phenotypic variation. However, in corroboration with Senthilvel et al. [58] and Wei et al. [76], QTL mapping studies on DH and RIL populations revealed that chromosome 3 holds a promising trait and was associated with nitrogen use efficiency (NUE) in rice. In P deficiency, several QTLs have been reported on chromosome 3 [55,60,61], on chromosome 5 [48,59,61], on chromosome 9 [59][60][61], and on chromosome 11 [48,60,61], which explained PV ranging from 4.4% to 16.6% in different mapping populations of BILs, ILs, RILs, and BC 2 F 3 in rice. Under low K, Wu et al. [39] identified a total 21 QTLs, related to PH, TN, SDW, RDW, relative potassium concentration in plant (RKC), relative potassium use efficiency (RKUE), and relative potassium uptake (RKUP) on chromosomes 2, 3, 7, and 8, that collectively explained PV from about 8% to about 15% using DH populations. On chromosome 3, qFGN_3.6 significantly associated with RTN and other traits, such as PH, SDW, RDW, and RKUP, were flanked by RG179-RG403 and RZ284-RZ394, contributing PV of 9.7% to 14.4% [39]. Thus, it is suggested that the identified QTL harbor regions support the reported NPK QTLs, and are refined with chromosomal regions of genes and QTLs that were involved in multiple stress tolerance and transportation of nutrient elements from soil to shoots and grain, which could be valuable information for the understanding of their genetic and physiological mechanisms in future molecular breeding programs in rice. Seeds of 230 ILs, parents, and four checks (PSB Rc82, NSIC Rc222, Apo, and IR74371-70-1-1) were sown in a seedling nursery bed, and 21-day-old seedlings were transplanted with a single seedling per hill. The early generation of the backcross population (BC 1 F 2 ) was grown in one generation under low-input, rainfed, and irrigation conditions during the 2011 wet season (WS), and was followed by four consecutive generations over the DS and WS during 2012 and 2013. Phenotypic screening for efficient selection was practiced based on higher grain yield under six nutrient conditions. Subsequently, this led to the development of 230 ILs. The detailed of the breeding strategies are described in Jewel et al [142]. The NuUE phenotyping experiment was laid out in an alpha lattice design with two replications, using a plot size of two rows × 12 plants/row, and a spacing distance of 0.2 × 0.2 m. NPK nutrients in the form of urea, superphosphate, and muriate of potash were applied at 160, 50, and 50 kg ha −1 in the DS, and at 90, 30, and 30 kg ha −1 in the WS, respectively. The checks and parents were replicated in all six nutrient conditions. The NPK fertilizers were applied five times in splits, as in basal level, and at 20, 36, 54, and 72 days after transplanting (DAT), respectively, against -N, -P, -NP, -NPK, 75N, and NPK conditions (Figure 4). Field management, including pest control, weeding, and irrigation, followed IRRI's standard experimental farm practices to avoid adverse effects on grain quality.

Measurements of Agronomic and Yield-Attributed Traits
Each IL was selected in the middle row of plants, with three replications, and the plants were randomly sampled from each plot for phenotypic evaluation of six agronomic and yield-related traits. The six respective nutrient conditions, -N (PK fertilizer), -P (NK fertilizer), -NP (K fertilizer), -NPK (zero fertilizer), 75N (75% NPK fertilizer), and NPK (all nutrients), were recorded for seven prominent NuUE traits: agronomic efficiency (AE), partial factor productivity (PFP), grain yield (GY), biomass yield (BY), filled grains per plant (FGN), 1000-grain weight (1000-Gwt), and percentage of spikelet fertility (PSPF). Data analysis of the key phenotypic traits was performed with R software [143], and SPSS 17.0 software (IBM, Armonk, NY, USA) was used for the analysis of descriptive means, analysis of variance (ANOVA), Pearson correlation coefficient, and regression analysis, which were calculated among the ILs.

Calculation of Agronomic Efficiency (AE) and Partial Factor Productivity (PFP)
The increase in grain yield for each kg of fertilizer applied is known as AE. Using the formulas below, we calculated AE and PFP along with nitrogen fertilizer use condition. AE (kg kg -1 ) = grain yield (N fertilized-N unfertilized) in kg ha −1 /Fertilizer N in kg ha −1 [144]. AE was computed for all the nutrient use efficient ILs, under six different nutrient conditions. The recommended fertilizer application rate was NPK = 160-50-50 kg ha −1 . Among them, three possible AE formulas were depicted: Similarly, PFP of applied nitrogen is also called nitrogen use efficiency. PFP (N) = YN (crop yield with applied N (kg ha −1 )/amount of fertilizer N applied (kg ha −1 ) [144]. Based on applied NPK fertilizer, the performance of the selected and fixed nutrient use efficient ILs, across six different levels of nutrient conditions, was determined by the PFP of applied nitrogen using the following formula:

Measurements of Agronomic and Yield-Attributed Traits
Each IL was selected in the middle row of plants, with three replications, and the plants were randomly sampled from each plot for phenotypic evaluation of six agronomic and yield-related traits. The six respective nutrient conditions, -N (PK fertilizer), -P (NK fertilizer), -NP (K fertilizer), -NPK (zero fertilizer), 75N (75% NPK fertilizer), and NPK (all nutrients), were recorded for seven prominent NuUE traits: agronomic efficiency (AE), partial factor productivity (PFP), grain yield (GY), biomass yield (BY), filled grains per plant (FGN), 1000-grain weight (1000-Gwt), and percentage of spikelet fertility (PSPF). Data analysis of the key phenotypic traits was performed with R software [143], and SPSS 17.0 software (IBM, Armonk, NY, USA) was used for the analysis of descriptive means, analysis of variance (ANOVA), Pearson correlation coefficient, and regression analysis, which were calculated among the ILs.

Calculation of Agronomic Efficiency (AE) and Partial Factor Productivity (PFP)
The increase in grain yield for each kg of fertilizer applied is known as AE. Using the formulas below, we calculated AE and PFP along with nitrogen fertilizer use condition. AE (kg kg −1 ) = grain yield (N fertilized-N unfertilized) in kg ha −1 /Fertilizer N in kg ha −1 [144]. AE was computed for all the nutrient use efficient ILs, under six different nutrient conditions. The recommended fertilizer application rate was NPK = 160-50-50 kg ha −1 . Among them, three possible AE formulas were depicted: Similarly, PFP of applied nitrogen is also called nitrogen use efficiency. PFP (N) = YN (crop yield with applied N (kg ha −1 )/amount of fertilizer N applied (kg ha −1 ) [144]. Based on applied NPK fertilizer, the performance of the selected and fixed nutrient use efficient ILs, across six different levels of nutrient conditions, was determined by the PFP of applied nitrogen using the following formula: (iv) PFP (N) = Y (+NPK) /F N (v) PFP (N) = Y (−P) /F N (vi) PFP (N) = Y( 75N) /F N Y +NPK = crop yield with applied NPK fertilizer (kg ha −1 ); F NPK = amount of fertilizer NPK applied (kg ha −1 ); F N = fertilizer N in kg ha −1 ; AE = agronomic efficiency applied nitrogen; PFP = partial factor productivity applied nitrogen.

Genotyping via 6K SNP Array
A total of 230 ILs of plant genomic DNA were extracted from leaf tissue using the cetyl trimethylammonium bromide (CTAB) method [145] and quantified by a NanoDrop 8000 spectrophotometer (ThermoFisher Scientific, Waltham, MA, USA). The concentration of DNA was adjusted to 50 ng µL −1 , and approximately 200 ng of DNA from each genotype used in the SNP array. The complete set of the BC 1 F 5 mapping population, along with the parents used for the genotyping with 6K SNP array technology, was followed by DNA quantification, incubation, and hybridization of bead chips, staining, and image scanning according to the manufacturer's instructions for the Illumina Infinium assay, and this work was conducted in the Genotyping Services Laboratory of the International Rice Research Institute. The resulting intensity data were processed using the genotyping module V2011.1 of Genome Studio software (Illumina Inc., San Diego, CA, USA) for SNP calling. The generated genotypic data indicated that 704 high-quality SNP markers were identified and further used in the construction of high-density linkage maps for NuUE QTLs, related to agronomic and yield-attributed traits.

Mapping of QTLs and Hotspot Regions for NuUE
The mapping of QTLs was carried out with IciMapping software (QTL IciMapping version 4.0) [146] using a single marker analysis (SMA) method. The permutation method was used to obtain an empirical threshold for claiming QTLs based on 1000 runs of randomly shuffling the trait values [147], and the logarithm of odds (LOD) value threshold for claiming QTLs was 2.5. The genetic distance (cM) between SNP marker positions was changed to physical distance (kb), with 1 cM equal to 260 kbp [148,149]. For seven promising critical agronomic and yield-related traits under six NuUE conditions, significant markers on the same chromosome, linkage disequilibrium (LD), and phenotypic trait data were calculated by the genetics package in R software [150]. A graphical representation of the linkage map was constructed using MapChart software [151]. Further, for fine-tuning of hotspot QTLs, we used the Q-TARO database (http://qtaro.abr.affrc.go.jp/) for the analysis of previously published co-localization of QTLs and genes related to soil stress tolerance in the rice genome.

Conclusions
In many plant breeding programs, the development of rice varieties with NuUE is being considered as a means to reduce the usage of NPK fertilizers by improving their use efficiency, thereby improving grain yield productivity. In the present study, selective introgression lines derived from Weed Tolerant Rice 1, as the recipient parent, and Hao-an-nong, as the donor parent, enabled us to identify a large number of QTLs (261 putative QTLs) linked with NuUE traits, of which 49 QTLs showed high PVE ranging from 20.25% to 34.68%. Among them, 22 QTLs reported as novel QTLs were responsible for PFP and AE traits. They were identified among the promising top four hotspot QTLs (QTLs harbor I-IV), which comprised more than ten QTLs associated with critical NuUE traits, such as 1000-Gwt, PFP, BY, FGN, GY, and PSPF. The hotspot regions of QTLs, expressed across all six NuUE conditions, suggested an underlying uniform basis of genetic mechanisms, contributing to the tolerance of these traits and associated with tightly linked genes or QTLs, or pleiotropic regulations. Corroboration with several research efforts, of the earlier discovered QTLs for NUE traits, revealed that most of the genomic region was confirmed and their positions precisely confirmed by significant SNP peak markers with high LOD and PVE values, and this could have potential value for introgression of the target using MAS. The list of QTLs and refined hotspot regions will facilitate further validation in systematic breeding for specific adaptability under low-input conditions, and suggest that hotspot genomic regions could be used as targets for a superior understanding of the NuUE mechanism and for improving NuUE traits in rice. In the future, identified promising harbor QTLs would be useful for developing elite introgression breeding lines comprising positive QTLs via marker-assisted selection, and also using them for carrying them forward by a pyramiding approach over the ILs, with desirable QTLs within the same population. Further, this will lead to fine-mapping and molecular cloning of the critical loci that will be useful for enhancing grain yield and quality under low-input fertilizer management conditions.