Identification of quantitative trait loci associated with nitrogen use efficiency in winter wheat

Maintaining winter wheat (Triticum aestivum L.) productivity with more efficient nitrogen (N) management will enable growers to increase profitability and reduce the negative environmental impacts associated with nitrogen loss. Wheat breeders would therefore benefit greatly from the identification and application of genetic markers associated with nitrogen use efficiency (NUE). To investigate the genetics underlying N response, two bi-parental mapping populations were developed and grown in four site-seasons under low and high N rates. The populations were derived from a cross between previously identified high NUE parents (VA05W-151 and VA09W-52) and a shared common low NUE parent, ‘Yorktown.’ The Yorktown × VA05W-151 population was comprised of 136 recombinant inbred lines while the Yorktown × VA09W-52 population was comprised of 138 doubled haploids. Phenotypic data was collected on parental lines and their progeny for 11 N-related traits and genotypes were sequenced using a genotyping-by-sequencing platform to detect more than 3,100 high quality single nucleotide polymorphisms in each population. A total of 130 quantitative trait loci (QTL) were detected on 20 chromosomes, six of which were associated with NUE and N-related traits in multiple testing environments. Two of the six QTL for NUE were associated with known photoperiod (Ppd-D1 on chromosome 2D) and disease resistance (FHB-4A) genes, two were reported in previous investigations, and one QTL, QNue.151-1D, was novel. The NUE QTL on 1D, 6A, 7A, and 7D had LOD scores ranging from 2.63 to 8.33 and explained up to 18.1% of the phenotypic variation. The QTL identified in this study have potential for marker-assisted breeding for NUE traits in soft red winter wheat.

Introduction Wheat (Triticum aestivum L.) is among the most widely grown crops in the world and accounts for roughly 20% of the global dietary calories and protein consumed by humans per annum [1]. It is therefore crucial to continually improve wheat productivity and quality to meet grain demands in an era of increasing human population. Since the onset of the Green Revolution in the mid-20th century, wheat yield gains in the eastern United States' soft red winter wheat growing region have been largely attributed to active selection for performance under intensified nitrogen (N) management conditions [2]. However, a majority of the N applied to agricultural crops in this region is not harvested and thus subject to loss from the plant-soil system [3]. The unharvested N is associated with emissions of nitrous oxides [4], runoff and leaching of nitrate [5], and the degradation of aquatic and terrestrial ecosystems [6][7]. Therefore, the development of wheat cultivars that more efficiently take up and utilize applied N provides a means of promoting grower profitability and ensuring environmentally sustainable increases in wheat production.
Genetic improvement of N use efficiency (defined as the ratio of harvested grain per unit N applied; NUE) [8] in wheat has been previously reported under a range of N conditions in Europe [9], Central America [10], and the United States [11]. While the extent of genetic improvement under differing N conditions varied by study, the authors are generally in agreement that direct selection under multiple N rates will accelerate genetic gains in NUE. This hypothesis was further supported by the detection of significant genotype by N rate (G × N) interactions for grain yield and N traits in recent studies of genotypic variation in wheat [9, [12][13]. However, the authors also note that direct selection under multiple N rates may be cost-prohibitive for wheat breeders. In response to this dilemma, Cormier, et al. [14] proposed identifying genomic regions associated with N response, known as quantitative trait loci (QTL), to enable more efficient cultivar selection. Through this approach, breeders can efficiently screen germplasm for genetic markers associated with N response to assist in the development of high NUE cultivars.
Previous studies of cereal crops have searched for novel NUE traits and alleles in adapted breeding materials [15], landraces [16][17], and wheat wild relatives [18]. While these authors have successfully identified QTL, genes, and genotypes conferring high NUE, additional sources of genetic variation likely still exist within currently unexplored germplasm. This hypothesis has withstood testing in European began the process of dissecting the genetic variation underlying NUE in United States wheat germplasm by conducting genome-wide association studies for N-traits in the central and eastern Unites States, respectively. It is therefore crucial to follow up these investigations using biparental mapping populations to validate findings and to enrich potentially rare NUE alleles for subsequent QTL analysis.
Successful QTL mapping for complex traits including NUE is dependent on the selection of suitable parents, evaluation of appropriate population sizes, multi-environment testing, and the development of high-density genetic maps The present study sought to: i) construct a high-density genetic map using GBS derived markers; and ii) identify and validate QTL associated with N-related traits under normal and reduced N conditions. We aimed to detect stable and impactful marker-trait associations for NUE and N-related traits with direct application for marker-assisted breeding programs.

Plant materials
The present study employed a population of 136 RILs, derived from a cross between 'Yorktown' and VA05W-151 (herein known as "YT×151"), and a population of 138 DHs, derived from a cross between parental lines Yorktown and VA09W-52 (herein known as "YT×52"

Experimental design
Both populations were grown under rainfed conditions in four Eastern Virginia testing environments (defined as a "site-season") described in Table 1. Testing sites included the Eastern Virginia Agricultural Research and Extension Center near Warsaw, VA in the 2015-2016, 2016-2017, and 2017-2018 winter wheat growing seasons (WR; 37˚99' N, 76˚78' W) and a commercial production field near New Kent, VA in 2017-2018 (NK; 37˚54' N, 76˚89' W). Seeds were treated with Raxil MD (triazole, Bayer Crop Science) and Gaucho XT (imidacloprid, Bayer Crop Science) in all testing environments to control diseases and insects, respectively. Foliar pesticides and herbicides were applied throughout the growing season in all environments to further mitigate pest pressure (S1 Table). Experimental units, seven-row yield plots that measured 2.74 m long × 1.52 m wide in WR and 4.88 m long × 1.78 m wide in NK, were sown at a seeding rate of 480 seeds m -2 to ensure consistent stand establishment. The two populations utilized a type-2 modified augmented design [43][44] in all environments due to limited seed and land availability. The design was comprised of 12 statistical blocks with each block consisting of 30 plots (three rows and 10 columns per block) to facilitate spatial adjustments in all testing environments. Each block consisted of a centrally positioned primary check line (OH08-161-78), two randomized secondary check lines (OH08-172-42 and 'Sisson', [PI 617053]) [45], and 12 randomized experimental lines with each check and experimental line being grown under two N rates in adjacent plots. Parental lines were replicated three times per population under each N rate per environment. Wheat lines in each statistical block were grown under high (134 kg N ha -1 ; HN) and reduced (67 kg N ha -1 ; LN) spring N rates that were foliar applied as liquid urea ammonium nitrate split over Zadoks growth stages 25 and 30 [46]. Pre-plant N, Phosphorous, potassium, and sulfur were applied according to soil test recommendations prior to planting (S1 Table).

Phenotypic measurements
In the WR site, anthesis date was recorded when half of the anthers had emerged from the florets. Maturity date was recorded in 17WR (where "17" refers to the growing season) and 18WR where it was defined as the date at which when 75% of the peduncles within a plot turned yellow. In all testing environments, plant height was measured before harvest from two random locations within the center rows of each yield plot. A 1.0 m above-ground biomass (AGBM) sample was cut from the center row of each plot at harvest maturity (grain moisture concentration � 160 g kg -1 ) in all testing environments. Biomass samples were oven dried at 60˚C for 72 hours and weighed to estimate AGBM yield for all experimental and check lines. Samples were then threshed to estimate harvest index (calculated as the ratio of grain biomass per unit of total AGBM) by separating grain from straw and chaff tissue and recording their respective weights. Grain and straw tissues were ground through a 2 mm sieve and homogenized to estimate tissue N concentration via combustion analysis using a Vario Max Cube elemental analyzer (Elementar Analysensysteme, Hanau, Germany). Results of the combustion analysis enabled calculation of grain N concentration, NUE, N uptake efficiency (NUpE), and N utilization efficiency (NUtE). Moll,et al. [8] defined the terms NUpE and NUtE as quotients of NUE where NUpE is calculated as the ratio of aboveground N at harvest per unit of N applied and NUtE as the amount of grain produced per total aboveground N in the plant at harvest. Plots were combine harvested (Wintersteiger Classic; Wintersteiger, Ried, Austria) in all testing environments at harvest maturity and adjusted to 0 g moisture kg -1 to determine grain yield.

Statistical analysis
An analysis of variance (ANOVA) was performed for the replicated parents using the lme4 package [47] in the R statistical computing environment [48]: Where the trait response (Y ijkl ) is a function of the overall mean (μ), the fixed effect of the ith wheat line (G i ), the fixed effect of the jth N rate (N j ), the random effect of the lth replication (R l ) nested within the kth environment (E k ), the interactions of the ith wheat line with the jth N rate and the kth environment (GN ij and GE ik ), the interaction of the jth N rate and the kth environment (NE jk ), their 3-way interaction (GNE ijk ), and the residual error (ε ijkl ). This was followed by an ANOVA for parental lines within environment as the effects of testing environment and its interactions were significant for most traits. Means comparisons were conducted using least significant differences for single effects and their interactions. Following the ANOVA, a bivariate correlation analysis calculated Pearson's correlation coefficients using trait means for progeny in each population under LN and HN conditions.

Genotyping
Genomic DNA was isolated from seedlings of RILs, DHs, and parents at the three-leaf stage using an LGC 218 Genomics Oktopure™ robotic extraction platform with sbeadex™ magnetic microparticle 219 reagent kits at the USDA-ARS Eastern Regional Small Grains Genotyping Center (Raleigh, NC, United States). Three technical replicates of each parent and single replicates for RILs and DHs were submitted for genotyping. Genotyping-by-sequencing was performed using an Illumina HiSeq 2500 following the protocol of DNA digestion with the restriction enzymes PstI and MspI [49]. Sequence reads were aligned to the 'Chinese Spring' v1.0 reference genome [33] using the Burrows-Wheeler Aligner v0.7.17-r1188 [50]. TAS-SEL-GBS v5 [51-52] was used to perform SNP calling, and to generate consensus calls for replicated parent samples. Resulting genotypic data was then filtered to retain SNPs with less than 20% missing data frequencies and less than 5% heterozygous calls using VCFTools v0.1.16 [53] and BCFTools v1.9 [54]. SNP calling and filtering was performed on each population separately. Missing data was not imputed and only SNPs with differing parental homozygous calls were retained. The remaining genomic DNA from each population was used to amplify polymorphic markers from a set of 116 SNP and simple sequence repeat (SSR) markers used in routine screening of the uniform and regional breeding nurseries (S2 Table).

Construction of genetic maps and detection of QTL
Polymorphic markers were used to construct high-density linkage maps for the YT×151 and YT×52 populations. The SNP calls were converted to an ABH parent-based format for the con-

Phenotypic variation and trait correlations
The ANOVAs for the replicated parental lines, N rates, and their interactions within each environment are shown in S3 and S4 Tables. Variance of parents and summary statistics of progeny for all traits at each N-environment for the YT×151 and YT×52 populations are reported in S5 and S6 Tables, respectively. Significant phenotypic variation was observed in the two populations under each N rate for a majority of the studied traits. All progeny reached anthesis over a seven and 10-day period in the YT×151 and YT×52 populations, respectively. The parents in both populations expressed similar plant heights, anthesis dates, and maturity dates within each testing environment despite variation at the Ppd-D1 locus in the YT×52 population (S2 Table). Parents of the YT×151 population differed significantly in NUE under LN supplies in all testing environments, while their NUE was not statistically different under the HN rate (Fig 1 and S5 Table). The YT×52 population expressed a similar trend but the parents only differed significantly under the LN supply in the 18NK environment (S6 Table).

Linkage map construction
After filtering for low-quality markers, linkage maps comprised of 3,918 markers spanning 2,962.8 cM in the YT×151 population and 3,147 markers spanning 2,491.7 cM in the YT×52 population ( Table 4). Within the YT×151 population, chromosome 6B had the lowest

QTL in the YT×151 population
The YT×151 linkage map was employed to detect QTL with a LOD score greater than 2.5 for 11 traits. A total of 66 QTL were identified for traits in one or more N-environments (S7 Table) of which 12 QTL were reproducible in two or more N-environments (   Table). Quantitative trait loci for nitrogen use efficiency in winter wheat

QTL in the YT×52 population
The YT×52 population linkage maps were employed to detect QTL in the same eight N-environments as the YT×151 population. Within this population, a total of 64 QTL were identified with a LOD score greater than 2.5 (S8 Table) and mapped to the A (25), B (19), and D (20) genomes. However, only eight of these QTL were reproducible in two or more N-environments; these were located on chromosomes 1A, 2D, 4A, and 7A (  (Fig 5) and did not co-localize with any non-reproducible QTL observed in this population (S8 Table).

QTL identified in both populations
Common single and reproducible QTL identified in both the YT×151 and YT×52 populations on 1B, 2D, 3B, 3D, 4A, 5D, 6A, and 7A (Tables 5 and 6 and S7 and S8 Tables). One common reproducible QTL identified in both populations was located within 10 cM of the FHB-4A locus on chromosome 4A. Additionally, the reproducible QTL found in the YT×151 population (QHi.151-3B) was located near a NUtE QTL that was identified on 3B in one N-environment (16WR-LN) in the YT×52 population. Similarly, the reproducible QTL for AGBM and NUE identified in the YT×151 population (QAgbm.151-6A and QNue.151-6A) were found in a similar position to an AGBM QTL from the YT×52 population in 16WR-LN on chromosome 6A. Six additional non-reproducible QTL that were common to both populations were identified on chromosomes 1B for NUtE (YT×151) and NUE (YT×52), 2D for grain N concentration (YT×151) and AGBM (YT×52), 3D for harvest index (YT×151) and maturity date (YT×52), 5D for grain N concentration (YT×151) and maturity date (YT×52) and on a second loci for grain N concentration (YT×151) and plant height (YT×52), and 7A for NUE in both populations.

Effects of QTL combinations on NUE
Reproducible QTL and their combinations were assessed for NUE under LN and HN conditions. Within the YT×151 population, the presence of QNue.151-1D and QNue.151-6A significantly increased NUE over testing environments under LN by 1.8 and 1.7 kg kg -1 , respectively  Table 7). The NUE QTL on 7D, QNue.151-7D, produced a significant increase in NUE under both LN (1.8 kg kg -1 ) and HN (1.2 kg kg -1 ) supplies. Nitrogen use efficiency further increased with two or three combined QTL from the YT×151 population. The YT×52 NUE QTL,  (Tables 7 and 8).  Quantitative trait loci for nitrogen use efficiency in winter wheat investigation of 225 European wheat lines, AGBM was further shown to be highly heritable (h 2 = 0.79) over environments [9], but traditional phenotyping proved to be exceptionally laborious. Frels, et al. [68] sought to overcome this constraint by identifying high-throughput vegetative indices that were predictive of AGBM to enable more efficient phenotyping. While the investigation showed some promise for selection of genotypes expressing high AGBM, the tested vegetative indices were not consistently predictive of AGBM over growing seasons. Future work is required to improve the predictive ability of current AGBM models but the trait appears to be a suitable target for NUE breeding in winter wheat.

Trait variation and associations in the mapping populations
The present study also identified a strong relationship between NUE and its component trait, NUpE, under LN and HN rates in both populations, while the other component trait, and 2) the influence of each trait is dependent on the environment, management practices, and genetic material being tested. It is therefore up to wheat breeders within specific regions to identify traits that most limit NUE and those for which phenotypic assessment is feasible.

Identification of QTL for NUE traits
As expected, the GBS genotyping platform increased the number of polymorphic SNPs available in the present study compared to previous investigations that utilized the 90K-SNP array for soft red winter wheat [73][74]. The use of GBS increases marker density through a whole genome scan approach as opposed to screening a select number of potential polymorphisms and makes results between investigations more reproducible due to the availability of both physical and genetic positions. The present study detected similar numbers of SNPs, chromosome coverage, and marker densities as a previous bi-parental population of winter wheat that utilized GBS [34]. However, map length was exceptionally low on chromosome 6B in the YT×151 population which may indicate low allelic diversity on 6B between the two parents. Previous wheat NUE mapping studies have identified QTL for N traits on every chromosome [22, 26, 75-77] but often lack reproducibility due to the quantitative nature of the studied traits. The low number of reproducible N-related QTL in these studies was fairly consistent with our ability to detect a large number of reproducible QTL within a single population in the present study. However, the number of reproducible QTL identified in the present investigation nearly doubles when QTL detected in a single N-environment are assessed over both populations. It therefore becomes challenging to identify which of these QTL have potential application for marker-assisted breeding programs that often employ fewer than 100 genetic markers. In response to this challenge, Quraishi, et al. [78] conducted a meta-analysis of QTL for NUE to identify 11 major chromosomal regions linked to N use efficiency. The NUE QTL, QNue.151-6A, identified in the present study was located near the QTL on 6A described in the aforementioned study that co-localized with a known glutamine synthetase gene (GS1) [65] and may therefore have a practical application in marker-assisted breeding due to its involvement in the assimilation of ammonium into amino acids. Additional QTL within proximity to the NUE locus on 6A were associated with NUtE and kernel weight per spike [79] and kernel weight [80] in previous studies of N response. The authors further linked this QTL on 6A to TaGW2, which influences grain size and weight [81][82]. The IWA4036 marker located near the reproducible QTL for NUE on chromosome 6A flanks a known FHB resistance locus [60]. However, the other flanking marker, IWA3483, was not segregating in the population and may therefore indicate that the genetic regions on 6A governing FHB resistance and improved NUE and N traits are different in the YT×151 population. Interestingly, this QTL was only significant under LN conditions in the YT×151 population and was also associated with AGBM under LN conditions in the YT×52 population and may therefore be a valuable marker for yield improvements under N limiting conditions. The QTL, QNue.151-1D, inherited from Yorktown was shown to increase NUE in multiple N-environments and AGBM in one N-environment in the VA05W-151 background. Bordes,et al. [83] similarly identified a QTL associated with NUE in bread wheat with close proximity to the NUE QTL, QNue.151-1D, identified in the present investigation. The high-molecularweight glutenin subunit gene, Glu-D1, was located more than 25 cM from this QTL and, therefore, is not likely the candidate gene [84]. However, Sun, et al.
[26] reported a QTL on 1D associated with root and shoot weight in wheat seedlings grown in hydroponic culture at a similar mapping interval as the QNue.151-1D and, therefore, may indicate a role in seedling vigor. QNue.52-7A was found within the same mapping interval as QTL previously linked to grain yield and spikes per m -2 [21] and a QTL associated with NUtE and spikes per m -2 [79] under low and high N conditions. A majority of the QTL identified in the present study validated findings from previous studies, yet QNue.151-7D was not identified in marker-trait association studies for N response and may represent a novel QTL. While previous mapping studies have identified QTL within proximity to QNue.151-7D, they are primarily associated with a vernalization allele, VRN3, located more than 50 million nucleotides from the presently identified QTL [19, 85, 86].

QTL co-segregating with known genes
In addition to evaluation for parental variation in N response, the parents shared similar alleles for vernalization (Vrn-A1, Vrn-B1, and Vrn-D1), photoperiod (Ppd-A1, Ppd-B1, and Ppd-D1), and plant height (Rht-B1 and Rht-D1). However, the Vrn-A1 and Ppd-D1 loci differed between Yorktown and VA09W-52 and the resulting allelic variation on 2D resulted in a strong effect QTL that explained relatively large percentages of the variation for NUE, grain yield, plant height, and anthesis date in multiple environments and N supplies. Indeed, this result is consistent with previous NUE investigations of wheat grown in the eastern [ The negative effects on NUE conferred by the QTL identified near FHB-4A was much higher than expected upon development of the mapping populations. Indeed, associations between FHB resistance and yield per se is not uncommon as previous investigations found a small negative impact on grain yield [90] and grain quality [91] conferred through a FHB resistance locus on chromosome 5A. The known yield penalties associated with FHB resistance have discouraged some winter wheat breeders from utilizing exotic sources of resistance including that of spring wheat variety 'Sumai-3' [92 -94]. However, the mechanisms governing a reduction in NUE by the QTL linked to FHB-4A in the present study requires further investigation to determine if this linkage can be broken.

Implications for breeding
The QTL identified in this study that were not linked to the FHB-4A or Ppd-D1 may merit use in marker-assisted-breeding programs as they were not associated with major physiological traits or known sources of disease resistance. Genetic markers for high NUE were previously reported near QNue.151-1D, QNue.151-6A, and QNue.52-7A and, therefore, offer the most potential for immediate marker deployment. The other NUE QTL, QNue.151-7D, requires validation in unrelated populations to determine its value in NUE breeding programs. All QTL identified in this study will also benefit from introgression into diverse backgrounds to gauge their overall value to soft red winter wheat breeding programs.
Supporting information S1