Genetic dissection of developmental responses of agro-morphological traits under different doses of nutrient fertilizers using high-density SNP markers

The production and productivity of rice (Oryza sativa L.) are primarily influenced by the application of the critical nutrients nitrogen (N), phosphorus (P), and potassium (K). However, excessive application of these fertilizers is detrimental to the environment and increases the cost of production. Hence, there is a need to develop varieties that simultaneously increase yields under both optimal and suboptimal rates of fertilizer application by maximizing nutrient use efficiency (NuUE). To unravel the hidden genetic variation and understand the molecular and physiological mechanisms of NuUE, three different mapping populations (MPs; BC1F5) derived from three donors (Haoannong, Cheng-Hui 448, and Zhong 413) and recipient Weed Tolerant Rice 1 were developed. A total of three favorable agronomic traits (FATs) were considered as the measure of NuUE. Analysis of variance and descriptive statistics indicated the existence of genetic variation for NuUE and quantitative inheritance of FATs. The genotypic data from single-nucleotide polymorphism (SNP) markers from Tunable Genotyping-By-Sequencing (tGBS) and phenotypic values were used for locating the genomic regions conferring NuUE. A total of 19 quantitative trait loci (QTLs) were detected, out of which 11 QTLs were putative on eight chromosomes, which individually explained 17.02% to 34.85% of the phenotypic variation. Notably, qLC-II_1 and qLC-II_11 detected at zero fertilizer application showed higher performance for LC under zero percentage of NPK fertilizer. The remarkable findings of the present study are that the detected QTLs were associated in building tolerance to low/no nutrient application and six candidate genes on chromosomes 2 and 5 within these putative QTLs were found associated with low nutrient tolerance and related to several physiological and metabolic pathways involved in abiotic stress tolerance. The identified superior introgressed lines (ILs) and trait-associated genetic regions can be effectively used in marker-assisted selection (MAS) for NuUE breeding programs.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Introduction significant yield advantage of 1.51 t ha -1 was achieved by reducing the N fertilizer application by 22.41% [20]. Hence, an organized and convenient procedure is required to establish a stage-specific fertilizer dose application to improve the nutrient balance for increasing agronomic advantages with the least environmental hazards. Three key favorable agronomic traits (FATs), plant height (PH), tiller number (TN), and leaf chlorophyll content (LC), and other yield-attributing traits contribute significantly to better plant architecture [36][37][38][39]. These traits were highly influenced by the application of a combination of complete N, P, and K fertilizers. Interestingly, these traits are controlled by the combined effects of several genes (polygenes) with a large influence of environment. Hence, the dissection of these traits that follow non-Mendelian segregation had become difficult and complex, especially when applying traditional breeding tools and methodologies. Recent advances in molecular marker and genomics technology paved the way for dissecting the genetic basis of the complex inheritance of key agronomic traits. Genetic analysis of quantitative trait loci (QTLs) can estimate and identify the influences of different genes responsible for quantitative traits with high statistical power. Further, the identified genes conferring these quantitative traits can be introgressed through precise MAS either singly or in combination with other useful traits/genes through markerassisted pyramiding approaches [21,40,41].
DNA-based molecular markers such as restriction fragment length polymorphism, amplified fragment length polymorphism, cleaved amplified polymorphic sequence, and simple sequence repeats were used to generate genotypic information to determine the genetic characteristics of the lines. However, the generation of information is time-consuming and laborintensive [42][43][44] as these markers are slower in work to provide accurate genotypic information for predicting candidate gene loci due to the low resolution of QTL mapping, environment-dependent expression of QTLs, and gene epistasis [44,45]. Recently developed highthroughput sequencing technologies are robust, accurate, and more informative for generating millions of single nucleotide polymorphism (SNP) markers in a shorter period [43,46,47]. These SNPs are most promising to develop genetic linkage maps for QTL dissection of many traits with higher chromosomal coverage [48][49][50][51]. To date, several QTLs have been reported for N, P, and K deficiency tolerance traits by using different QTL mapping methods in different genetic backgrounds of MPs such as recombinant inbred lines (RILs), backcross inbred lines (BILs), introgression lines (ILs), doubled haploids (DHs), chromosome segment substitution lines (CSSLs), and BC 2 F 3 families in rice [10,41,52]. According to a comprehensive literature survey and exploring information available in the Gramene database (http://archive. gramene.org), more than 150 QTLs for N, 130 QTLs for P, and 15 QTLs for K have been associated with deficiency tolerance traits in rice. These significant QTLs are located on eight chromosomes (1, 3, 4, 5 7, 8, 9, and 12) for N [21,[53][54][55][56][57][58][59]. Similarly, P deficiency tolerance QTLs have been reported on five chromosomes (1, 2, 6, 11, and 12) [10,41,53,[60][61][62][63][64][65][66][67]. K deficiency tolerance QTLs are located on three chromosomes (3, 5, and 8) [68,69]. Notably, rather than a deficiency of the individual nutrient, either N or P or K, the combination of the complete N, P, and K deficiency has a significant impact on agronomic traits such as tillering ability, PH, and LC, which are key contributors for determining total grain yield [70,71].
To the best of our knowledge, the identification of QTLs conferring individual N, P, and K deficiency tolerance has been less explored. Interestingly, there was no comprehensive information on QTL data regarding the effect of fertilizer application in different growth interval (DGI) stages and during complete deficiency of N, P, and K fertilizers in rice. Hence, understanding the genetic information of N, P, and K deficiency tolerance of FATs under DGI stages of fertilizer doses is crucial to finding suitable rice cultivars for N, P, and K use efficiency (NuUE). In the current study, we used a set of early backcross MPs derived from a cross of three donors, Haoannong (HAN), Cheng-Hui 448 (CH448), and Zhong 413 (Z413), with a recipient parent, Weed Tolerant Rice 1 (WTR-1), to study the genetic variation for deficiency tolerance for N, P, and K at different stages. In addition, we detected genomic regions conferring N, P, and K defficiency tolerance by using phenotypic evaluation under three different doses of N, P, and K fertilizer, 338.90 kg ha -1 (100%), 271.10 kg ha -1 (80%), and 0 kg ha -1 (0%), and genotyping by tGBS technology to detect the QTLs associated with the trait of low input tolerance. The major objectives of the study were to (i) identify the FAT responses under three different doses of fertilizer application, 100%, 80%, and 0% for N, P, and K conditions; (ii) understand the genetic basis of FAT responses at DGI stages under three different doses of fertilizer application through correlation and heritability estimation; (iii) identification of promising tolerant ILs; and (iv) detect the main-effect QTLs (M-QTLs) at DGI stages under different doses of fertilizer by also estimating the additive effects of M-QTLs (v) enumeration of the candidate genes within the major effect QTLs based on the in silico candidate genes analysis using online databases.

Plant materials and developing mapping populations
The parental line WTR-1 developed by using a novel Green Super Rice breeding strategy, was used as the recurrent parent while HAN, CH448, and Z413 were used as donor parents [72][73][74] for the development of MPs for the identification of the genomic regions controlling N, P, and K use efficiency. A set of 120, 67, and 56 BC 1 F 5 plants were generated from the HAN × WTR-1 (DP-1), CH448 × WTR-1 (DP-4), and Z413 × WTR-1 (DP-7) cross combinations to make a total of 243 BC 1 F 5 plants. The parental lines and their MPs along with four checks (Rc 222, SL8, Rc 192, and PSB Rc 82) for each trait underwent phenotypic evaluation to assess the genetic variation and dissect the locus conferring the FATs at DGI stages of fertilizer application under three different N, P, and K fertilizer doses.

Experimental layout
The introgression lines and their parental lines were seeded on the experimental farm at IRRI, Los Baños, Philippines (14˚11N, 121˚15E), during the wet season of 2017. A total of six individual plots separated by ridges were modified in a grid to maintain two replications. An alpha lattice experimental design was followed to analyze the experimental significance of the laidout experiment. Each plot size was 245.76 m 2 and each entry was planted in 0.4-m rows, with 0.2-m distance between rows and space of 0.2 m between adjacent plants. The 400-μm-thick transparent polyethylene sheets were inserted into the soil to cover the individual plots. At 21 days after seeding, seedlings were transplanted in an area of 491.52 m 2 in two replications. The whole experiment was laid out with three different rates of N, P, and K fertilizer. The first rate of fertilizer consisted of complete N, P, and K that is also a recommended dose of fertilizer designated as 100% (N-P 2 O 5 -K 2 O = 100:30:30 kg ha -1 ); the second rate is a suboptimal dose, designated as 80% (N-P 2 O 5 -K 2 O = 80:30:30 kg ha -1 ); and the third rate is zero fertilizer dose, designated as 0% (N-P 2 O 5 -K 2 O = 0:0:0 kg ha -1 ). Furthermore, during crop growth, the fertilizers were applied in four DGI at 5, 20, 35, and 50 days after transplanting (DAT). Table 1 presents the detailed DGI at different fertilizer doses. stem tiller was counted as one tiller and included in the total tiller number, whereas LC was measured in the middle leaf portion by using the portable, non-destructive chlorophyll meter SPAD-502Plus (Konica Minolta, Japan). The mean LC value indicates the relative chlorophyll content. The phenotypic data were collected from three plants selected from the middle of each entry at four DGI stages every 15 days after fertilizer application. The first fertilizer application was done at 5 DAT and was designated as stage-I. Similarly, stage-II, stage-III, and stage-IV DGIs were designated for fertilizers applied at 21 DAT, 36 DAT, and 51 DAT. The phenotypic data were expressed as mean values for each of the observations recorded. The mean values were further used for the localization of genomic regions conferring these traits.

Genotyping
Genomic DNA was extracted from individual plants of the MPs along with the parental lines using the DNeasy Plant Mini Kit (Quiagen, USA). The quality and quantity of the isolated DNA were analyzed using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA) and sent to Data2Bio for tGBS analysis, in which all samples were sequenced through 10 Ion Proton runs. Further, SNP calling and filtration of the generated sequencing information were processed by adopting the methodology developed by Data2Bio, LLC (https://www.data2bio.com), as described by Ali et al. [73]. For the SNP filtration, the threshold of 50% missing rate (LMD50 = minimum calling rate) was set to filter out SNPs having more than 50% missing rate across the MPs. The missing SNPs that were showing � threshold missing rate were imputed using Beagle software ver. 3.3.2 [75]. Additionally, another level of SNP filtering was performed to remove mono and hetero SNPs for the individuals of the MPs using TASSEL ver. 5 and finally exported as HapMap format. The obtained HapMap format SNP marker file was again processed in MS-Excel ver. 2010 to prepare an input file for the QTL analysis [76]. The complete work-flow of tunable GBS data analysis and its use for QTL mapping studies are represented in Fig 1.

Construction of linkage map and QTL analysis
The inclusive composite interval mapping (ICIM-ADD) method is available in QTL IciMapping ver. 4.1 software [77] for scanning the FAT data at DGI stages from three cross combinations with different nutrient fertilizer doses (0%, 80%, and 100% N, P, and K fertilizer). Before construction of the linkage map, SNPs were analyzed for segregation distortion using the Chi-square test at a significance of P = 0.05. The SNPs that were showing significant segregation distortion and were co-located in the same physical position were filtered out using the redundant marker removal feature of the software. The processed high-quality SNP markers were used for constructing the genetic map to assign markers on each linkage group. Kosambi mapping function was used to compute the genetic distance using recombination fraction (cM) [78]. The precisely estimated phenotypic data of FATs at different DGIs for each nutrient dose (0%, 80%, and 100% NPK) for each cross combination and high-resolution SNP marker data were used for the QTL analysis. In order to identify the QTLs with high precision and to declare significant QTLs, a 1000 permutation test at the 0.01 level of significance was considered. The linkage map and QTL mapping were constructed for each MP derived from each cross combination (DP-1, DP-4, and DP-7) separately. The detected QTLs which showed >10.00% R 2 value were considered as major QTLs.

Candidate gene analysis
To identify the candidate genes within the QTLs, the interval region of both flanking marker positions was considered. The physical position of the SNP flanking markers was used to determine the size of each QTL in kb. Publically accessible MSU Rice Genome Annotation Project (Osa1) release 7 (http://rice.plantbiology.msu.edu/) and the Rice Annotation Project (https://rapdb.dna.affrc.go.jp/) were explored to select the candidate genes residing in each detected QTL. Further, in silico gene expression analysis carried out using the genome wide expression-profiling database, RiceXPro (http://ricexpro.dna.affrc.go.jp). This is the repository of the microarray profiling data of Nipponbare cultivar isolated from different growth and developmental stages.

Statistical analysis
The phenotypic data of each trait at four DGI stages for the three fertilizer doses were analyzed using analysis of variance (ANOVA) for assessing the significance of the experiments at the level of significance P = 0.05. DMRT and Fisher's t-test were used to ascertain the significant difference between the genotypes and to compare with the check values. The correlation between the traits at different fertilizer doses (0%, 80%, and 100% N, P, and K) was computed by Pearson's correlation analysis and the correlation matrix was visualized using corrplot package in R (https://github.com/taiyun/corrplot) [79]. The phenotypic mean values of all the traits were prepared using MS-Excel ver. 2010 and the analysis was carried out using PBTools (Version 1.4, http://bbi.irri.org/products). Chi-square goodness of fit was used to analyze the segregation pattern of SNP markers before using these for the QTL analysis. The multiple regression model of maximum likelihood was employed for the composite interval mapping.

The genetic variation of phenotypic traits among different MPs
The early backcross MPs derived from three different genetic backgrounds were evaluated for three FATs under 100% NPK (338.90 kg ha −1 ), 80% NPK (271.90 kg ha −1 ), and 0% NPK (0 kg ha −1 ) treatments at four DGI stages (I, II, III, and IV). The segregation pattern of FATs was determined based on the skewness and kurtosis values. For the MP DP-1, normal distribution was observed for the trait PH at stage II under 100% fertilizer rate; at stages I, II, and III at 80% fertilizer rate; and at stages I and II at 0% fertilizer rate. For the MPs DP-4 and DP-7, most of the traits at different stages and different fertilizer rates showed normal distribution (Table 2 and S1 Fig). ANOVA showed significant variation existing in the three concentrations of N, P, and K fertilizers at specific growth stages for each trait. The F-test statistic values showed the existence of significant variation at the level of significance P = 0.05, indicating a large amount of genetic variation among the MPs (  in MP DP-1. As expected, the lowest CV values were recorded in the second stage of LC in DP-1, DP-4, and DP-7 in all the different doses of fertilizer. The broad-sense heritability (H 2 b ) was estimated for each trait under different doses of NPK fertilizers at DGI stages. Under 100%, H 2 b ranged from 15% to 51%. The highest H 2 b was noted in TN-IV (51%), followed by LC-II and PH-IV, being the same estimated value of 46%. Similarly, in 80% and 0% NPK conditions, H 2 b ranged from 10% to 65% and 6% to 65%, respectively.

Pearson correlation matrix (PCM) analysis
The PCM was performed in DGI stages in three genetic backgrounds of MPs to assess the cor-

SNP discovery and construction of saturated linkage map
All the genotypes, including parental lines, underwent genotyping by using the advanced and cost-effective genotyping platform tGBS. It was used for the extraction of SNP markers for the localization of the genetic factors influencing FATs under different fertilizer doses. The Table 3.

Analysis of variance (ANOVA) and its components for all the FATs in DGI stages under different doses of fertilizer.
Fertilizer dose Components PH-I  TN-I  LC-I  PH-II  TN-II  LC-II  PH-III  TN-III  LC-III  PH-IV  TN-IV  LC-  genotyping platform tGBS yielded extensive sequence information from sequencing three MPs. A total of 10,432, 14,117, and 7,865 raw reads were obtained from the MPs as DP-1, DP-4, and DP-7, respectively (Table 4). A novel SNP calling and filtering method was followed to eliminate missing values (�50 missing data, LMD 50) from the tunable GBS pipeline (Data2-Bio, LLC) (https://www.data2bio.com), as reported in Ali et al. [73]. In the initial SNP filtering, as many as 4,666, 5,957, and 3,147 SNPs were isolated from MPs DP-1, DP-4, and DP-7, respectively, whereas, in the next level of filtering, 4,286 SNPs from DP-1, 3,622 SNPs from DP-4, and 2,610 SNPs from DP-7 were isolated (Table 4). Notably, out of these SNP markers, 35.10% in DP-1, 44.64% in DP-4, and 47.67% in DP-7 were redundant and excluded before further analysis. Finally, the SNPs were filtered out at the final level of filtering based on the segregation distortion regions in the genetic linkage map. A total of 57.80% of the SNPs demonstrated genetic distortion (P = 0.05) in the DP-1 population, followed by 44.63% in DP-4 and 38.72% in DP-7. The SNP markers obtained from all the filtration levels were used for linkage map construction and linkage analysis. In total, a 953.71 cM chromosomal spanning region was covered by 2,782 SNPs in DP-1, 1,398.42 cM by 2,005 SNPs in DP-4, and 834.11 cM by 1,361 SNPs in DP-7. In each population, an average of 230, 167, and 113 SNP markers were distributed across all chromosomes. The highest number of polymorphic markers was present on chromosome 1 (324 SNPs) while the lowest was on chromosome 9 (172) in DP-1, whereas, in DP-4 and DP-7, the highest number of markers was located on chromosome 1 (279) and chromosome 2 (223). However, chromosome 9 (70) and chromosome 4 (11) had the lowest number of polymorphic markers available in DP-4 and DP-7, respectively (Table 5 and Fig 3). The average distance of each SNP marker varied from 0.13 to 0.41 Mb, whereas the longest distance recorded was 1.47 Mb on chromosome 4 in the population DP-7.

Molecular mapping of genomic regions conferring FATs
In order to locate the QTLs influencing the traits PH, TN, and LC at DGI stages under different fertilizer doses, QTL analysis was performed using the ICIM-ADD methodology. A total of 19 main-effect QTLs (M-QTLs) were identified using the composite interval mapping model of QTL IciMapping under LOD threshold of 3.00. The significant QTLs were identified on chromosomes 1, 9, 10, and 12 in 100% NPK conditions; on chromosomes 1, 2, 3, 5, and 12 in 80% NPK conditions; and on chromosomes 2, 3, 4, 5, 6, 8, and 11 in 0% NPK conditions. The detected M-QTLs explained phenotypic variation ranges from 1.89% to 34.85% with LOD score ranges from 3.02 to 28.70. The number of QTLs per trait ranged from one to four, and the highest numbers of QTLs (four QTLs) were associated with LC-II in populations DP-1, DP-4, and DP-7 (S1 Table). A total of four QTLs were discovered for LC-II, TN-I and II, and

Putative QTLs with high statistical power
To detect QTLs with higher statistical power, an extreme threshold parameter, 1000 permutation test at P = 0.01, was considered. A total of 11 putative QTLs (�15% PVE) were detected for the three FATs at DGI stages under three different fertilizer doses (100%, 80%, and 0% NPK). The QTL qPH-III_12 was the only one detected at 100% NPK for the FAT, PH at stage III on chromosome 12 between the markers S12_852867 and S12_2965504. Under suboptimal fertilizer dose, three QTLs for PH, qPH-III_2, qPH-IV_3 and qPH-I_5,were detected on chromosome 2, 3 and 5, while another single QTL for LC (qLC-II_1) and TN (qTN-III_5), were detected on chromosome 1 and 5. Simialrly in zero fertilizer conditions, five QTLs (qTN-III_3, qPH-I_4, qTN-I_5, qLC-II_8 and qLC-II_11) for PH, TN and LC were detected on five chromosomes 3, 4, 5, 8, and 11. However, the significant M-QTLs were identified only from MP DP-1, except four QTLs detected in the cross combination CH448 × WTR-1 (DP-4) and Z143 × WTR-1. The distribution of detected significant main-effect QTLs under three different NPK conditions is depicted in Fig 4. Additionally, minor effects QTLs were also detected   for all the traits. The list of total minor effect QTLs including the information of flanking markers along with their LOD score, PVE and the additive effect of parental source of the allele are presented in S1  -II_1), all of them were identified in a single MP, DP-1 (S1 Table).
Similarly, QTL analysis under zero fertilizer dose revealed that a total of five main-effect QTLs were detected for the traits of PH-I, TN-I and III, and LC-II. Three QTLs (qTN-III_3, qPH-I_4, and qTN-I_5) were detected on chromosomes 3, 4, and 5 and they explained PV of 17.02%, 19.17%, and 28.44% at stages I and III of fertilizer application. Another two QTLs (qLC-II_8 and qLC-II_11) were associated with a single trait (LC) at first stages of fertilizer application. QTLs on chromosome 3, 4, and 5 were contributed by the negative additive effect from donor alleles, whereas QTLs on chromosomes 8 and 11 were contributed by the positive additive effect from the WTR-1 allele.

Exploring putative QTLs for candidate genes for future breeding programs
The putative QTLs detected from the extreme threshold parameters were further explored to showcase the possible candidate genes within the loci to understand the molecular and physiological mechanisms underlying the traits conferring NPK deficiency tolerance. A total of 19 M-QTLs were located on all chromosomes, except on chromosome 7. Further, the possible candidate genes within the M-QTLs were filtered based on the SNP flanking marker positions with a threshold score of �1 Mb interval regions of M-QTLs ( Table 6). The lowest number (3) and highest number (25) of candidate genes were associated with a single trait (TN at stages III and IV) in the MP DP-1. A total of 120 genes were identified (S2 Table), of which 78.33% (94 genes) were functionally annotated while 21.66% (26 genes) were reported as hypothetical, retrotransposon proteins, and expressed proteins. The 94 functionally annotated genes were classified into six biological functions related to cellular component (CC), biological process (BP), and molecular function (MF) (Fig 5). Out of these 94 genes, 63 candidate genes were involved in MF, CC, and BP, followed by 18 genes for BP and MF, 9 genes for BP and CC, 3 genes for CC, and a single gene for BP.

Discussion
For the ever-increasing global population and to meet food demand, the development of rice varieties with higher grain yield is essential [80]. The application of the key nutrients N, P, and K either organically or through chemical fertilizer plays a foremost role in increasing yield and sustaining soil fertility [81]. With the recent trends in crop breeding, higher rice productivity has been successfully attained by applying high fertilizer doses [9,35,82,83]. However, without knowledge of the correct stage, timing, and dose of fertilizer application, any assurance for Genetic dissection of developmental responses of agro-morphological traits under different fertilizer doses increasing yield cannot be determined. The excessive use of fertilizer is a major contributor to increasing soil, water, and environmental pollution, along with farm operation costs [9,18,33,84]. China is leading in per hectare usage of fertilizer (300 to 350 kg ha -1 N in Jiangsu Province), and this amount is around four times higher than the average world fertilizer consumption for rice production [18,33]. Hence, it is crucial to undertake a systematic breeding program involving identifying genotypes with higher input use efficiency and genomic regions derived from these lines to improve elite varieties to assure sustainable crop production. This would not only reduce environmental and human hazards but would also improve the livelihood of farmers by reducing farm operational costs along with a higher expected grain yield. Hence, in the present study, we have developed MPs generated from three different donors (HAN, CH448, and Z413) and evaluated them under three different fertilizer doses, 100% (222.23 kg ha -1 N), 80% (177.77 kg ha -1 N), and 0% (0 kg ha -1 N), to study the genetic variation for low input use efficiency, especially for N, P, and K. Further, we identified the genomic regions influencing FATs at DGI stages under different fertilizer doses while dissecting the molecular genetic information.

Dissecting molecular genetics of nutrient deficiency tolerance in rice
ANOVA carried out using different MPs for FATs at different stages indicated the existence of significant variation among all the MPs for the target traits. As expected, all the genotypes, including parents of the MPs, showed a higher mean performance for all the AFTs at 100% NPK, followed by 80% and 0% NPK. This trend observed mainly because of the higher doses of the key nutrients N, P, and K, which are proven to enhance rice yield. Several researchers also reported the similar effect of key nutrients on plant growth and yield-attributed traits [85][86][87][88]. The variance components of FATs indicated that the segregation pattern at different DGI stages under different fertilizer doses varied from normal to skewed distribution ( Table 2). This clearly confirmed the influence of few or many genes with a cumulative and additive effect, which is difficult to dissect using traditional low-resolution genotyping platforms. Hence, to understand the molecular genetic basis of the traits that influence the genomic regions conferring FATs at different DGI stages under varied fertilizer doses, a high-resolution and informative genotyping platform was employed. Numerous high-quality SNPs retrieved from the advanced genotyping tool tGBS were used for the molecular mapping of the key traits. A total of 19 M-QTLs were identified from three different MPs for the three FATs at four DGI stages under three different fertilizer doses. Out of these QTLs, 13 QTLs from DP-1, 4 QTLs from DP-4, and 2 QTLs from DP-7 explained PV ranging from 1.89% to 28.98%, 7.74% to 28.44%, and 30.62% to 34.85%, respectively. With the extreme threshold parameters, eight M-QTLs were considered as putative QTLs possess �1 Mb QTL regions, whereas the remaining 11 QTLs possess �1 Mb genetic regions of QTLs ( Table 6). The putative QTLs qPH-III_2, qLC-I_2, qPH-I_5, qTN-III_5, qTN-IV_5, qLC-II_8, qTN-I_9, and qTN-III_12 suggested that these are closely associated with the respective traits (PH, TN, and LC) at DGI stages of fertilizer application. However, the majority of these QTLs on chromosomes 2, 5, and 8 were attributed to WTR-1 allele, whereas the QTLs on chromosomes 9 and 12 by HAN.

Potentiality of the putative QTLs and promising ILs
The mean performance of all the individuals of the MPs showed higher values than their parents for all the FATs (  [13,[89][90][91]. However, in the present study, we could identify several ILs that showed ) from DP-7, have been commonly identified as promising ILs that show constitutively increased LC at DGI stages. These ILs, being highly responsive to fertilizer application, could be an imperative source of NuUE traits for future breeding programs and also an excellent source for the genetic dissection of tolerance of low inputs in rice. These potential ILs are highly useful for medium to low input marginal farmers who may not incur much of their cost in buying fertilizer for rice cultivation. Therefore, the area and production of rice would increase, which could further help in ensuring food security and farmers' livelihood. This is one of the remarkable findings and applications of the present study and could be attributed to the contribution of the putative QTLs detected from the potential genetic backgrounds.

Promising QTLs and comparison with previous QTLs related to deficiency of N, P, and K fertilizers
The rice Gramene database (http://archive.gramene.org), QTL Genome Viewer (http://qtaro. abr.affrc.go.jp), and previous reports from comprehensive literature surveys showed hundreds of QTLs associated with morphological, physiological, and biochemical traits that influence individual nutrient fertilizer deficiency in rice [10,41,52,92]. Among the 19 M-QTLs detected in the present study, 14 M-QTLs were detected in the same genomic region that was previously reported. These QTLs were significantly associated with more than 20 traits reported under N, P, and K deficiency tolerance, and were located on all chromosomes except chromosomes 7 and 8. The remaining five QTLs were found to be novel (  N). Some of these QTLs were co-localized with low N, P, and K tolerance on chromosomes 2, 4, 6, 10, 11, and 12. These QTL clusters might play an important role in NuUE in rice, and their related SNP markers could be useful for MAS in low input breeding programs.
By comparing with the previously reported QTLs, most of the QTLs were sharing common genetic regions at DGI stages in all populations. The QTLs were associated with other multiple QTLs, ranging from one to 11 under N, P, and K deficiency tolerance (Table 7). In the zero percentage of fertilizer condition, a total of three putative QTLs related to relative PH (cm) (qRPH), relative shoot dry weight (g/p) (qRSDW), and relative potassium uptake (mg/p) (qRKUP) shared the same genomic region with qLC-I on chromosome 2. Among the three QTLs, qRPH had the highest phenotypic variation of 14.70%, with an LOD score of 4.04, identified by using 123 DH lines under low-K stress conditions [68]. On chromosome 3, six QTLs related to N concentration in leaf sheaths plus stems (qNS%_3.1), N concentration in leaf blades (qNL%_3.2), dry weight of roots (qDWR_3.3), grain density per panicle (qGD_3a), number of filled grains per panicle (qFGP_3a), and number of spikelets per panicle (qSP_3b) shared the same genetic region with qTN-III_3 at 29.07-28.88 Mb [58,93]. These QTLs were responsible for nitrogen use efficiency (NUE) and were also associated with yield-attributed traits under low nitrogen rates. Two genetic regions on chromosome 4 at 15.43-21.82 Mb and 31.28-32.42 Mb region of two QTLs (qPH-I_4 and qTH-II_4) were shared with three QTLs related to shoot and root growth traits under low nitrogen and phosphorus tolerance and another two of them were related to phosphorus transporters [57,67]. By using the F 9-10 generation of 169 recombinant inbred lines (RILs) derived from a cross between IR64 (O. sativa L. var. indica), and Azucena (O. sativa L. var. japonica), Thi et al. [93] identified 44 QTLs for 15 Genetic dissection of developmental responses of agro-morphological traits under different fertilizer doses agronomic and NUE-related traits on all chromosomes, except chromosome 9. Among these traits, seven NUE-related traits such as agronomic NUE, number of tillers, total fresh matter, dry weight of roots, total dry matter, dry weight of sheaths plus stems, and dry weight of leaf blades were shared with the current study for two QTLs (qTN-I_5 and qTN-IV_5) on chromosome 5 at 17.80-17.86 Mb region and 17.80-18 Mb region, respectively. However, the highest number of QTLs (11 QTLs) for low nitrogen and phosphorus tolerance were shared with qTN-II_6 at the 2.37-10.28 Mb region on chromosome 6. Earlier, under low nitrogen conditions, Hu et al. [94] identified three QTLs for nitrogen content in shoots by using 116 DH populations, which have been developed through anther culture of F 1 hybrids from indica rice variety Taichung Native 1 (TN1) and japonica rice variety ChunJiang 06 (CJ06). One of the main-effect QTLs, qNCm6-12 (nitrogen content in plant shoots at mature stage), is significantly associated with the present QTL related to TN at stage II of fertilizer application. These QTLs mapped between RM527 and RM3 and explained PV of 9.73%. In addition, Feng et al.
[58] identified 28 QTLs for yield-attributed traits on seven chromosomes under low N conditions. Two QTLs (qPL-6a and qSP-6a) for panicle length and number of spikelets per panicle explained PV of 15.58% and 6.40% from the analysis of 138 F 14 RILs, which were derived from a super hybrid rice (Xieyou 9308) in China. Another QTL related to absorbed NUE (qaNUE) was also located in the same position [93]. The remaining seven QTLs were associated with low P tolerance QTLs. Of these, qPUP_6 was associated with P uptake [95], four QTLs (qRRDW, qRSDW, qRTDW, and qRPUC) for root traits and P uptake in rice [64], and two QTLs (qRRL, and qRWRSR) for root growth and weight traits [67]. Interestingly, to date, there are no reported QTLs matched to/shared with the currently identified major QTL qLC-II_8 (PVE of 22.85%) at the 8.21-8.92 Mb region on chromosome 8, which indicates that these QTLs are novel loci controlling for LC. However, Shimizu et al. [62], Wang et al. [96], and Tong et al. [97] reported 11 QTLs associated with LC on chromosomes 1, 2, 3, 4, 7, and 12 under low N and P conditions. As compared to the above studies, the current study identified four QTLs for LC on chromosomes 8 and 11 that are novel genetic regions under zero fertilizer conditions. The individual QTLs for LC had PVE ranging from 10.79% to 30.62%.The alleles from WTR-1 were in the direction of increasing the LC. These results indicate that these DGI stages of fertilizer application and cluster QTLs were significantly associated with low N, P, and K tolerances in rice.
Under suboptimal doses of fertilizer, a total of six QTLs (qLC-II_1, qPH-III_2, qPH-IV_3, qPH-I_5, qTN-III_5, and qTN-III_12) were associated with LC, TN, and PH located on five different chromosomes (1, 2, 3, 5, and 12), explaining an average PV of 19.09%. The smallest genetic region of qPH-III_2 (0.11 Mb) was located on chromosome 2, and qPH-IV_3 (1 Mb) was located on chromosome 3. Two QTLs, qPH-I_5 (1.04 Mb) and qTN-III_5 (0.33 Mb), were located on chromosome 5, and those are not shared with any reported QTLs. Therefore, these QTLs were considered as novel loci for DGI stages of fertilizer application. The remaining two QTLs, qLC-II_1 on chromosome 1 and qTN-III_12 on chromosome 12, at 6.15-24.88 Mb region and 0.85-2.97 Mb region, were shared. A total of three QTLs were related to low N and two QTLs to low P tolerance traits. Three QTLs (qTNCS_1b, qSY_1b, and qGD_1a) under low N on chromosome 1 shared the same genetic region with LC-II_1. The first two QTLs for total nitrogen content of shoot and straw yield were noticed by Cho et al.
[54] by using 166 F 8 lines derived from a cross between a Korean tongil type of rice, variety Dasanbyeo, and a Chinese japonica variety, TR22183. Similarly, Feng et al. [58] identified three QTLs for grain density (GD) per panicle by using 138 F 14 RILs, under low N. One of the QTLs, qGD_1a (6.04% PVE) located on chromosome 1, was shared with qLC-II_1. Two other QTLs (qRTDW_1 and qRSDW_1) were shared with qLC-II_1 on chromosome 1. These two QTLs were associated with low P tolerance. Lastly, chromosome 12 at the 0. 85-2.97 Mb region showed the presence of two QTLs (qRTHK_12_1 and qRDW_12) for root thickness [100] and root dry weight [99] by using RILs under low N conditions. On the other hand, for the recommended dose of fertilizer, four QTLs (qLC-II_1, qTN-I_9, qTN-II_10 and qPH-III_12) located on four different chromosomes (1, 9, 10, and 12) had considerable association with earlier reports of several QTLs related to low P and N tolerance in rice.
For yield-attributed and NUE traits related to QTLs under low N conditions, Lian et al. [99] and Dai et al. [102] reported a total of 15 main-effect QTLs for root and shoot weight and 13 QTLs distributed on all 12 chromosomes. Of these, two QTLs, for relative root weight, qRRW_1a [99], and seminal root length, SRL_1.2 [101], on chromosome 1 shared a common genetic region with qLC-II_1. Similarly, a single QTL, qDRW-9.5 [100], on chromosome 9 was shared with qTN-I_9 under low N conditions. The genetic region of 4. 70-14.56 Mb associated with qTN-II_10 is shared with two QTLs: for nitrogen accumulation amount, qNAA_10 [102], and phosphorus uptake, qPUP_10 [53]. Lastly, chromosome 12 at the 5.91-6.59 Mb region is associated with qPH-III_12 (PVE of 23.06%), which is not shared with any previous reports of QTLs and is considered to be a novel locus. In summary, among all the shared genetic regions, QTLs detected on chromosomes 1, 4, 5, 10 and 11 confer N and P deficiency tolerance, those on chromosomes 3 and 4 confer N deficiency tolerance, those on chromosome 4 confer P deficiency tolerance, and those on chromosome 2 confer K deficiency tolerance in rice. The clusters of QTLs detected under different doses of fertilizer at DGI stages and novel QTLs are valuable genetic resources to identify useful genes and introgression in elite lines to improve rice in NuUE.

Putative candidate genes and functions
The development of highly nutrient-responsive genotypes for low nutrient application is imperative during this era when the area and production of rice are under threat. MAS is one of the main strategies for accelerating the process of developing lines with higher NuUE. Effective MAS using major-effect QTLs mainly depends on the highest co-segregation of the markers and trait. Understanding of the molecular mechanisms and gene structure of the loci detected increases the efficiency of MAS. Hence, in the present study, we explored the candidate genes within the major-effect loci interval region (�1 Mb) and their predicted molecular and physiological functions. A total of eight putative QTLs detected on five different chromosomes (2, 5, 8, 9, and 12) were explored for the candidate gene analysis survey (Table 7). A total of 120 genes were identified in the promising QTL interval regions. Each QTL has a variable number of underlying associated genes ranging from 3 genes (58.2 kb) on chromosome 3 to 30 genes (711.90 kb) on chromosome 8. In silico analysis revealed that 78.33% of the genes are functionally annotated while the remaining genes are expressed and hypothetical proteins. The list of all possible candidate loci and their functions is detailed in S2 Table and Fig 5. These genes were functionally related to numerous physiological and molecular functions such as photosynthetic rate, synthesis of chlorophyll precursors, phosphate transporters (peptide transporter), zinc transporters, growth stimulators of auxin-responsive genes, abscisic acid (ABA) signaling pathways, and activation of various transcription factors, all of which were found within the promising M-QTL interval regions.
A QTL (qPH-II_2) located at 28.88-29.07 Mb on chromosome 2 associated with PH at stage II of fertilizer application was neighboring the candidate gene Os02g46140, which encodes for F-box protein domain involved in several biological functions, is related to phytohormone signaling pathways, and regulates various developmental processes as a part of the abiotic stress response mechanism [103]. Another candidate gene, Os02g46090, is associated with calcium signaling pathways. Several families of calcium sensors have been reported and they have played major roles in nitrogen metabolism and abiotic and biotic stress response mechanisms [104]. Recently, researchers identified a novel mutant of the calcium-dependent protein-kinase-encoding gene, esl4, which is significantly expressed in roots, shoots, and enzyme activity of several genes related to nitrogen metabolism [105]. Interestingly, two candidate genes, a tonoplast-localized low-affinity nitrate transporter, OsNPF7.2 [106], and a transcriptional regulator, GROWTH-REGULATOR FACTOR 4 (OsGRF4), along with a growth inhibitor of DELLA proteins [107], are located in the same genetic region on chromosome 2. The hotspot region of the root-specific transporter plays a vital role in intracellular allocation of nitrate in roots, especially in sclerenchyma, cortex, and stele, and ultimately leads to influences on plant growth. The interaction of GRF4 and DELLA proteins is involved in physiological activities to regulate multiple genes for carbon and nitrogen metabolism, plays a significant role in the homeostatic coordination of nitrogen metabolism, and also increases shoot biomass and PH [107,108]. Overexpression and mutant analysis of OsNPF7.2 showed that the significant expression of OsNPF7.2 caused an increase in TN by regulating the cytokinin and strigolactone pathway [109]. Another QTL, qLC-I_2 at the 30.45-30.56 Mb region, was found to be located on chromosome 2, adjacent to another candidate gene, OsGS1;1 (cytosolic glutamine synthatase), which regulates nitrogen metabolic pathways and also influences plant growth and development [110,111]. The GS1 proteins have a major role in generating glutamine, which is the primary form of remobilized nitrogen during natural senescence in leaves for long-distance transport [112]. In addition to that, qLC-I_2 is adjacent to the diacylglycerolacetyl transferase (OsDGAT) (Os02g48350) gene and interacts with OsTCP19 in regulating seedling establishment by modulating stress signaling molecules and abscisic acid pathways in various abiotic stresses [113].
Three QTLs, qPH-I, qTN-III, and qTN-IV (17.47-18.00 Mb), were found to be located on chromosome 5 adjacent to Os05g30970, and they encode for copine-like protein. It is mainly involved in the nutritional role of glutamate for aiding in seedling establishment under nitrogen deficiency [114]. In the same position are three other possible putative candidate loci: Os05g30870 (OsRLCK185), for diverse roles in plant growth and development and stress responses [115]; and Os05g30240, encoding for pentatricopeptide-repeat protein for upregulation in salt stress conditions [116]. Hence, to overcome nutrient deficiency, balanced mineral nutrients are essential for optimal plant growth and development. Each of these genes/loci played a critical role in various physiological and molecular responses such as cytokinins negatively regulating Pi (phosphate) starvation and also regulating metabolic changes under low N [117,118], and phytohormones such as auxins and seconday metabolites that are involved in the maintenance of homeostasis, root hair development, and signaling pathways under low NPK tolerance mechansim [119,120]. Similarly, phytohormones, ABA, jasmonates (JA), and their associated biosynthetic genes at vegetative and reproductive stages regulate various signaling pathways, leading to adaptation to nutrient deficiency through the root architecture, and maintaining N, P, and K homeostasis [121][122][123]. The in silico expression analysis was performed for the candidate genes, Os02g46140, Os02g46090, and Os02g48350 on chromosome 2 and Os05g30970, Os05g30870, and Os05g30240 on chromosome 5 located within M-QTL regions using RiceXpro [124]. The expression analysis indicated wide differential expression pattern of these candidate genes (S2 Fig). The high level of expression was observed in ovary development for the two candidate genes, Os02g46090 and Os05g30240, whereas for the other genes high level of expression was noticed in spikelet hull (lemma and palea) during the early stages of seed development and higher expression during roots and panicles development. Hence, these candidate genes within the M-QTL interval region can be considered as the promising putative candidate genes. However, further validation and developing functional markers are necessary for their effective application in breeding programs.

Conclusions
Genetic dissection of low N, P, and K tolerance at DGI stages of fertilizer application is an important target area for modern breeding programs in rice. Hence, in the present study, we identified the genomic regions that confer low N, P, and K tolerance for favorable agronomic traits at different stages of fertilizer application and selected ILs possessing the highest tolerance of low NPK using three different MPs derived from three donors (HAN, CH448, and Z413). ANOVA and descriptive statistics indicated the existence of significant genetic variation among the FATs and the predominance of non-Mendelian inheritance. To identify the genetic factors that influence these quantitatively inherited traits, QTL analysis was carried out by using the precisely estimated phenotypic values and high-quality SNPs derived from the tGBS genotyping platform. A total of 19 M-QTLs on different chromosomes were detected. Among these M-QTLs, eight QTLs were considered as putative QTLs, with the smallest locus size (�1 Mb) contributing the highest toward trait expression. Among the putative QTLs, qTN-I_92, qTN-III_5, qLC-I_2, and qLC-II_8 detected under 100%, 80%, and 0% NPK were found to be responsible for transgressive segregation of the ILs for the trait. Notably, qLC-I_2, qTN-IV_5, and qLC-II_8 detected at zero fertilizer application showed higher performance for LC under 0% of NPK fertilizer. These QTLs not only help in building a tolerance of low N, P, and K nutrient simultaneously but also improve genotypes to make them highly responsive to lower nutrient application. This is one of the remarkable achievements of the current study, which helps low-input and marginal farmers to cultivate rice without incurring high fertilizer cost and to eventually ensure food security and sustainable agriculture. On the other hand, in silico functional annotation of candidate genes within the putative QTLs indicated that two and five candidate genes found to confer tolerance of low N, P, and K and related to several physiological and metabolic pathways were also found to be involved in abiotic stress tolerance. However, additional investigation is needed for further confirmation to examine the potential physiological and molecular mechanisms. These studies can help in understanding the underlying complex genetic interactions involved in nutrient use efficiency and the identification of more efficient breeding materials containing these genetic factors. Furthermore, the detected genomic regions related to stage-specific tolerance of low fertilizer doses and promising ILs can be useful for MAS and future breeding programs for low-input conditions. Supporting information S1