QTL Mapping and Candidate Gene Identifying for N, P, and K Use Efficiency at the Maturity Stages in Wheat

Nitrogen (N), phosphorus (P), and potassium (K) are the three most important mineral nutrients for crop growth and development. We previously constructed a genetic map of unigenes (UG-Map) based on their physical positions using a RIL population derived from the cross of “TN18 × LM6” (TL-RILs). In this study, a total of 18 traits related to mineral use efficiency (MUE) of N/P/K were investigated under three growing seasons using TL-RILs. A total of 54 stable QTLs were detected, distributed across 19 chromosomes except for 3A and 5B. There were 50 QTLs associated with only one trait, and the other four QTLs were associated with two traits. A total of 73 candidate genes for stable QTLs were identified. Of these, 50 candidate genes were annotated in Chinese Spring (CS) RefSeq v1.1. The average number of candidate genes per QTL was 1.35, with 45 QTLs containing only one candidate gene and nine QTLs containing two or more candidate genes. The candidate gene TraesCS6D02G132100 (TaPTR gene) for QGnc-6D-3306 belongs to the NPF (NRT1/PTR) gene family. We speculate that the TaPTR gene should regulate the GNC trait.


Introduction
Wheat (Triticum aestivum L.) is a major crop around the world. It is a staple food for 30% of the world's population [1] and accounts for approximately one-fifth of the total calories consumed by humans [2]. Nitrogen (N), phosphorus (P), and potassium (K) are the most essential and important nutrient elements for plant growth and development [3]. N fertilizer is the main input used to boost grain yield and grain protein content, and agricultural production relies on synthetic N fertilizer to increase productivity [4]. P is an important component of plant biomolecules such as DNA, RNA, ADP, ATP, and NDPH. This makes it the second-most essential nutrient after N for plant development and crop production. As an essential macronutrient, K is involved in most biochemical reactions in plants, such as stomatal regulation, osmoregulation, protein synthesis, carbohydrate metabolism, and enzyme activation. In these processes, K can regulate the rate of generative reactions [5]. In recent decades, global fertilizer application has resulted in a significantly increased application of N, P, and K [6]. In an economic sense, fertilizer application has exceeded the optimal application rate, resulting in a massive waste of resources and economic losses [7]. Reliance on chemical fertilizers is clearly not a long-term solution, as both their production and application cause pollution and threaten environmental sustainability.
Mineral use efficiency (MUE) consists of two components: mineral uptake efficiency (MUpE), the ability of plants to remove minerals from the soil, and mineral utilization efficiency (MUtE), the ability of plants to use minerals to produce grain [8][9][10]. In wheat, the MUE of N, P, and K varies greatly among genotypes by growth stage, and the analysis of the MUE genetic basis is of great value [11][12][13][14][15][16]. The development of new high-yielding crop varieties under relatively low N/P/K environments is an effective way to improve MUE [17][18][19][20].
The traits associated with MUE are quantitative traits. Quantitative trait loci (QTL) analysis provides an efficient way to breakdown complex traits into loci and study their relative effects on specific traits [21]. In wheat, a large number of QTLs related to MUE have been mapped in low and high nutrient environments [22][23][24][25][26][27][28]. For example, Zhang et al. [29] detected 121 and 130 QTLs for N use efficiency traits under hydroponic and field treatments, respectively. Yang et al. [30] detected twelve new QTLs influencing root and biomassrelated traits for N uptake under three hydroponic N levels. Su et al. [31] identified seven QTLs related to P utilization efficiency under high and low P treatments. Yuan et al. [14] detected 55 and 68 QTLs related to P use efficiency under hydroponic and field treatments, respectively. Gong et al. [32] detected a total of 127 QTLs related to K use efficiency under 10 different N and P concentrations. Safdar et al. [33] identified two QTLs related to P use efficiency in a genome-wide association study and found that six genes were highly expressed in pyruvate metabolism and the TCA cycle under low phosphorus treatment.
We previously constructed a set of recombinant inbred lines (RILs) derived from the "TN18 × LM6" cross. The objective of the present study was to perform QTL analysis for MUE traits in field trials of RILs at different treatments of N, P, and K and then to find relatively stable QTLs and candidate genes.

Plant Materials
The population used for QTL analysis in this study consisted of a set of RILs derived from the cross of "TN18 × LM6" (TL-RILs, 184 lines) [14,29]. The female parent, TN18, is a cultivar developed by our group that was approved by the Crop Variety Approval Committee (CVAC) of Shandong Province in 2008. TN18 was planted on approximately 300,000 hectares per year in the Huang-Huai Winter Wheat Region, China. The male parent, LM6, is an elite wheat line whose mother is a sister line of the famous cultivar Jimai 22. The two parents showed obvious differences in MUE, with the LM6 having a higher MUE than the TN18 as a whole. A total of 184 lines of RILs were randomly selected from the original 305 lines to be studied.

Experimental Design
We constructed eight nutrient plots of 110 m 2 (10 × 11 m) in the experimental station of Shandong Agricultural University for mineral nutrient element experiments. The plots were divided by 1.5-meter-deep cement brick walls. The soil structure was maintained as a natural field with fertile soil [8]. Mineral nutrients were depleted by growing wheat and maize each year until nutrient levels met trial requirements. Three field trials were conducted during the 2012-2013, 2013-2014, and 2014-2015 growing seasons.
Four treatments were designed in the three trials: CK (target grain yield: 9000 kg/ha) and low N (LN), low P (LP), and low K (LK) (target grain yield: 6000 kg/ha for each corresponding element) [14,29]. All P 2 O 5 and K 2 O and 60% of the N were applied before sowing, and 40% of the N was applied at the stem elongation stage. Each trial was a complete block design with two replications in two nutrient plots. Each line of TL-RILs was sown in two rows, with 1 m in length and 25 cm between rows. Twenty seeds for one row were sown in 2012-2013 and 2013-2014, and 40 seeds for one row were sown in 2014-2015. Seeds were sown on 15 October, and plants were harvested on 13-15 June.

Trait Measurement
A total of 18 MUE traits were investigated (Table 1)

Data Analysis
Analysis of variance (ANOVA) was calculated using SPSS 18.0 software (SPSS Inc., Chicago, IL, USA). Since the traits were measured with mixed samples, a two-factor ANOVA was used. The genotypes and treatments were considered two factors using the average values of three growing seasons. All factors involved were considered sources of random effects. The generalized heritability (h B 2 ) was estimated by the following formula: , where σ g 2 is the genotype variance and σ e 2 is the total error variance.

QTL Detection and Candidate Gene Identification
The genetic map of unigenes (UG-Map) for the TL-RILs was used for QTL analysis. Each line of the TL-RILs was sequenced using RNA-Seq technology. The clean reads were mapped to the Chinese Spring (CS) RefSeq v1.1. A total of 31,445 polymorphic subunigenes were obtained, which represented 27,983 unigenes. The sub-unigenes were used to construct the UG-Map based on their physical positions [34]. Three software programs were used for QTL analysis: Windows QTL Cartographer 2.5 (http://statgen.ncsu.edu/qtlcart/ WQTLCart.htm (accessed on 1 August 2019), IciMapping 4.1 (http://www.isbreeding.net/ accessed on 28 January 2019), and MapQTL 6.0 (https://www.kyazma.nl/index.php/mc. MapQTL/ accessed on 1 August 2019) [35]. IciMapping 4.1 adopts the complete interval mapping method. Both the Windows QTL Cartographer 2.5 and MapQTL 6.0 software used the composite interval mapping method. For the three software programs, the LOD value was 3.0 and the step size was 0.5 cM. For a single software program analyzing data from the same treatment, we defined a stable QTL when it was detected over the AV + 1 environments. Finally, we employed the stable QTLs that were detected in at least two of the three software programs.
The unigene(s) covered or rounded (not exceeding the middle distance between two unigenes) by the peak position region of a QTL were regarded as the candidate gene(s) of the corresponding QTLs.

Phenotypic Variation and Heritability
For the parents of the RIL population, TN18 and LM6 showed obvious differences in most investigated traits (Table S1). For the RIL population, a wide range of variations existed. Transgressive segregation was observed for all of the 144 trait-treatment combinations (including AV). All of the 18 investigated traits in each treatment showed continuous distributions, which indicates their quantitative nature. The h B 2 for the investigated traits ranged from 61.75% (StKC) to 88.33% (GNC) ( Table S1).

Candidate Genes for Stable QTLs
A total of 73 candidate genes (including ncRNAs) for the 54 stable QTLs were identified (Table 2, Figure 1) that were covered or rounded by the regions of the peak positions of the corresponding stable QTLs. Among these candidate genes, 50 were annotated in Chinese Spring (CS) RefSeq v1.1 [2], eight were annotated in TL-RILs, and 15 were ncRNAs. The average number of candidate genes per QTL was 1.35 (73/54) with 45 QTLs (83.3%) containing only one candidate gene, seven QTLs (13.0%) containing two candidate genes, and two QTLs (3.7%) containing 3-8 candidate genes.
For the N/P/K accumulation traits, 11 candidate genes were identified, with four genes annotated in RefSeq v1.1, three genes annotated in TL-RILs, and four ncRNAs (Table 2, Figure 1). Of these, only two HC genes were obtained: TraesCS5D02G295200 for StNA and TraesCS6A02G372600 for StPA.

Discussion
The QTL locations were carried out mainly using the genetic map of DNA markers, and a great number of QTLs for MUE traits of N/P/K were identified [14,22,[25][26][27][28][29][30][31][32][33]. Using a genetic map of high-density DNA markers, a QTL could cover dozens of genes [36]. Therefore, it is difficult to identify highly reliable candidate genes for corresponding gene cloning. In this study, we used UG-Map for the RNA sequencing of the RIL population [34], which can directly determine the candidate genes for the QTLs. We identified 73 candidate genes from 54 QTLs of MUE traits, with an average of 1.35 candidate genes per QTL. This should provide convenience for the cloning of related genes and the genetic improvement of wheat breeding programs. In addition, the QTLs and their candidate genes in this study were identified according to the physical positions of the reference genome. It is difficult to accurately compare QTLs obtained based on molecular markers.
Because the cloned genes and their functions for MUE were not reported, the relationship between our candidate genes and MUE needs further verification. An example with phenotypic function is the candidate gene TraesCS6D02G132100 from QGnc-6D-3306 for GNC. The TraesCS6D02G132100 gene is annotated as a peptide transporter (TaPTR gene). NRT1 and PTR belong to the same gene family NPF (NRT1/PTR, nitrate transporter 1/peptide transporter) in terms of sequence homology, i.e., proteins that mediate the transmembrane transport of substances such as nitrate and small molecular peptides of 2-3 amino acids [37]. Therefore, we speculate that the TraesCS6D02G132100 gene should regulate the GNC trait.
TraesCS6A02G372600 was the candidate gene for QStpa-6A-10015 for StPA. This gene is annotated as a phosphate regulon sensor protein (TaPhoR gene). The PhoR protein can act as a protease in the inducible high-affinity phosphate transport system (Pst), which is composed of Pho regulators that accomplish phosphate transport behavior via transmembrane signaling [38]. Therefore, we speculate that the TraesCS6A02G372600 gene may regulate the StPA trait in wheat.
The TraesCS6B02G012600 gene was the candidate gene for QGpc-6B-361 for GPC. This gene is annotated as an F-box family protein (TaPFB gene). In Arabidopsis, FBX2 is a protein containing WD40 and F-box structural domains that can be used to recognize specific ubiquitination targets. For example, under phosphorus-sufficient conditions, the fbx2 mutant has higher levels of phosphoenolpyruvate carbolicase kinase 1 (PPCK1) transcripts and significantly higher root hair numbers and total intracellular phosphorus content than the wild type [39]. Therefore, we speculate that the TraesCS2D02G124400 gene may associate with GPC in wheat.
Using the QTL location of UG-Map for TL-RILs, we previously cloned two genes, the XIP (Xylanase Inhibitor Protein) gene, which regulates the SDS-sedimentation value and dough stability time [40], and the TaDHL (ATP-dependent DNA helicase) gene, which regulates the plant height [41]. These results further showed that the candidate genes in this study are favorable for corresponding gene cloning.
For TL-RILs, the male parent TN18 is derived from the cross "Laizhou137 × Yannong 19" and the female parent LM6 is derived from the cross "924402 × Lvmai13". The strain 924,402 is one parent of the cultivated variety "Jimai22". Of these parents, Yannong19 and Jimai22 are the authorized varieties with great area and the core parents in wheat breeding programs in China. Moreover, the parent of TN18 is also a core parent in China. Hundreds of authorized varieties derive from Yannong19, Jimai22, and TN18. Therefore, the QTLs and their candidate genes in this study may be widely identified in the varieties of China.