Introduction

Bread wheat (Triticum aestivum L.) is a primary staple cereal crop for human consumption with a global harvest of over 732 million MT from over 224 million hectares in 2015 (USDA-ERS 2016). Approximately 68% of the world wheat crop is used directly for food, and 21% is fed to livestock (Food Outlook: Global Market Analysis 2012). The U.S. Census Bureau projects that world population will grow from the present 7.3 billion to 9.4 billion by 2040 (United States Census Bureau 2016), which will increase the demand for staple food crops. Gilland (2006) estimates that the present per capita global average cereal production can be maintained through 2050, but to do so will require a 50% increase in the use of nitrogen (N) fertilizer. Cereal crops with improved nitrogen use efficiency (NUE) could increase grain yield at a given level of N fertilizer.

Globally, 12.7 MMT of N were applied as fertilizer to wheat, while 38.1 MMT were applied to all crops (FAO 2006). Therefore, on a global basis, wheat accounts for nearly one-third of crop fertilizer use. From 1990 to 2010, average fertilizer use per acre of wheat in the United States was 75 kg ha−1, and total average national use on wheat was 1.8 million nutrient tons (Economic Research Service 2013). Fertilizer cost has increased dramatically since 2000, at an annual rate of nearly 10% each year. U.S. farm price for N fertilizer in 2012 tripled since 2002. Not surprisingly, U.S. fertilizer application rates on wheat have not risen since 1988. Consequently, wheat genotypes that efficiently capture and convert soil N into grain protein are key to sustainably meeting rising global demand for grain and grain protein.

Nitrogen use efficiency (NUE) is defined classically as the yield of grain per unit of available N (Moll et al. 1982). NUE is the net result of N capture (uptake efficiency) and N conversion (utilization efficiency) (Moll et al. 1982; Sylvester-Bradley and Kindred 2009). To assess N capture and conversion, four component measurements are required: grain yield (GY), grain N concentration (GNC), straw yield (SY, above-ground, non-grain biomass), and straw N concentration (SNC). The components of grain protein yield are GY and grain protein concentration (GPC), and these two traits are the two principle targets of most hard wheat improvement programs. Both parameters respond to N fertilization, yet it has been difficult to improve both yield and protein simultaneously. The negative relationship between grain yield and GPC is well established and has a genetic basis (Oury and Godin 2007; Bogard et al. 2010; Gaju et al. 2011; Oury et al. 2003; Barraclough et al. 2010; Simmonds 1995). Among 27 wheat genotypes grown in 27 environments ranging in available N (Bogard et al. 2010), mean grain yield variability accounted for 78% of mean GPC variability across all environments with a slope of −0.015% g−1 m2. The two principle hypotheses to explain the negative correlation of grain yield and GPC are competition for energy resources by C and N assimilation processes (Munier-Jolain and Salon 2005) and a N dilution effect by carbohydrates (Acreche and Slafer 2009).

Grain N can be partitioned as coming from two sources: remobilization of N from vegetative tissues (N remobilization, NREM) and N assimilated by the plant after anthesis (post-anthesis N uptake, PANU). In winter wheat 15N tracer studies, 71.2% of grain N originated from remobilization of pre-anthesis N (Kichey et al. 2007). Similarly, using subtraction analysis, Bogard (2010) determined that 76–80% of winter wheat grain N originated from remobilization, and genetic variability for NREM was low. Post-anthesis N uptake and NREM are distinct processes. Although both PANU and NREM contribute to NUE, the traits were strongly negatively correlated (r = − 0.51 to −0.87) (Bogard et al. 2010). Grain protein deviation (GPD) is the deviation of observed GPC from the concentration predicted by the grain yield-GPC relationship. PANU, but not NREM, was significantly correlated (r = 0.67) with GPD. Selection for positive GPD has been proposed as a strategy for simultaneously improving grain yield and protein in wheat breeding (Monaghan et al. 2001; Guttieri et al. 2015).

Elite European winter wheat germplasm has demonstrated a trend toward improved nitrogen utilization efficiency (NUtE; grain yield per unit available N) and nitrogen harvest index (NHI) over the last 25 years (Cormier et al. 2013). The levels of applied N fertilizer (as high as 250 kg ha−1) and the crop yields (mean grain yields 5.8–9.0 t ha−1) in these trials, however, were much higher than typical for the Great Plains of the United States. Some of the purposes of this study, therefore, were to characterize the genetic variation for NUE traits within Great Plains hard winter wheat germplasm and to evaluate trends in the germplasm with year of release. Previous studies in European germplasm and growing environments also have demonstrated that semi-dwarf stature (Gooding et al. 2012; Loddo and Gooding 2012) and flowering time (Bogard et al. 2011b) significantly affect NUE traits. Therefore this study also explored the relationships of plant height and flowering time with NUE parameters.

Phenotyping for NUE traits is laborious and costly; therefore linked molecular markers would be advantageous for breeding. QTL analyses to identify genomic regions associated with NUE have been conducted using both biparental recombinant inbred and doubled haploid populations from Asian and European germplasm in both seedling and full-season experiments (An et al. 2006; Guo et al. 2012; Sun et al. 2013; Xu et al. 2014; Laperche et al. 2006, 2007, 2008; Zheng et al. 2010; Fontaine et al. 2009; Habash et al. 2007; Garcia-Suarez et al. 2010; Bogard et al. 2011b, 2013). Association mapping (Bordes et al. 2013; Cormier et al. 2014) approaches also have been applied to NUE traits. These QTL analyses were recently reviewed by Cormier et al. (2016) and consistently indicate a large number (32–380) of small-effect QTLs for NUE traits. Similar studies have not been reported in North American winter wheat germplasm. Therefore this study used genome-wide association analysis with high-density SNP markers to identify genomic regions that may be associated with NUE traits in Great Plains hard winter wheat germplasm.

Materials and methods

Germplasm

The test genotypes were a panel of 299 genotypes (258 hard red winter, 41 hard white winter) originating from both public and private wheat breeding programs in Colorado, Kansas, Montana, Nebraska, North Dakota, Oklahoma, South Dakota, and Texas. The panel was selected by breeders in the Great Plains region based on significance in production and/or contribution to the pedigrees of contemporary germplasm. Test genotypes in the panel included 193 cultivars and 106 breeding lines that were in advanced stages of testing in 2008. The test genotypes, their programs of origin, approximate release date, and pedigree are listed in Guttieri et al. (2015). Historical germplasm was included in the panel: the landrace, Turkey, the two ancestral cultivars, ‘Cheyenne’ and Kharkof, and five cultivars (‘Comanche’, ‘Wichita’, ‘Kiowa’, ‘Bison’, and ‘Tascosa’) developed before 1960.

Experimental design

Within each trial, nitrogen treatments (two rates) were applied as main plot treatments with two replications. Subplots of genotypes were arranged in 15 augmented incomplete blocks of two check genotypes (‘Jagger’ and ‘Settler CL’) with 20 test genotypes. These check genotypes also were included as test genotypes within the trial. The genotype Wichita was included twice in each trial. Therefore 60 plots of each of check genotype and 4 plots of each test genotype were grown in each environment.

Agronomic management

The experiment was grown at the University of Nebraska Agricultural Research and Development Center near Ithaca, NE. The soil type was a Sharpsburg silty clay loam, and soil test data were reported previously (Guttieri et al. 2015). Plots of each genotype were sown on Oct 4 2011 and Sept 25 2012 at approximately 260 seeds m−2 as four 3.2 m rows spaced 30 cm apart. To ensure adequate stand establishment, 1.8 cm of irrigation was applied to the 2013 trial on Oct 2 2012. Trials otherwise were rain-fed. Cumulative precipitation was reported previously (Guttieri et al. 2015). Trials were conducted under weed-free conditions, and in both years, multiple fungicide applications were made to control diseases. Details of pesticide applications are included in Online Resource 1.

Nitrogen rate treatments were established based on soil test data from samples collected in the spring of each year with the objective of providing contrasting levels of nitrogen availability (low, moderate) relative to a standard commercial target. Soil samples were collected to a depth of 25 cm on Mar 15 and 26 in 2012 and Apr 5 2013 in each of the 120 check plots in each year. Concentrations of macro- and micronutrients, pH, and organic matter in soil samples was determined by Ward Labs (Kearney, NE). Nitrogen application rate was established based on target yields of 4000 kg ha−1 in 2012 and 4700 kg ha−1 in 2013. Target yield was adjusted upward in 2013 based on improved stand establishment, higher soil pH and soil P concentration, and a long, cool spring with abundant rain. Using the South Dakota State University (SDSU) fertilizer recommendation for wheat (Gerwing and Gelderman 2005), which is based on soil test N for the top 60 cm of soil and projects a crop requirement of 0.0417 kg ha−1 N per kg ha−1 target yield, target N rates in the top 60 cm of soil were 167 and 196 kg ha−1 in 2012 and 2013, respectively. From the standard commercial calculations for the region (Ward 2016), residual nitrate N in the top 60 cm was estimated to be 70.8 and 93.8 kg ha−1 in late Spring of 2012 and 2013, respectively (Table 1). Fertilizer treatments were applied as 46% granular urea with a tractor-mounted drop spreader (Barber Engineering Company, Spokane, WA). In 2012, 15.6 kg ha−1 N was applied in the spring to the low N main plot and 60.5 kg ha−1 N was applied to the moderate N main plot on Mar 26 2012 when wheat was at early to mid-tillering (Zadoks stage 23–26). In 2013, no additional N was applied to the low N main plot, and 43.8 kg ha−1 was applied to the moderate N main plot on Apr 29 2013 when wheat was at mid-tillering. Residual + applied N fertility in low N main plots was 52 and 48% of target N rates in 2012 and 2013, respectively. Residual + applied N fertility in moderate N main plots was 79% and 70% of target N rates in 2012 and 2013, respectively.

Table 1 Spring residual and applied soil nitrogen and estimated total N available (NTAmax) to low and moderate N treatments in 2012 and 2013
Full size table

Phenotypic measurements

Anthesis or flowering date (FLO) was recorded visually as the date of anther exertion from 50% of main spikes. Physiological maturity date was recorded visually as the date of 90% chlorophyll depletion from the peduncle. Plant height (HT) was recorded at full stature as the height to the terminal spikelet less awns. Above ground biomass was collected at anthesis from 60 cm of row and 30 cm of row in 2012 and 2013, respectively. Samples were dried in forced air dryers at 55 °C for at least 7 days and weighed. At physiological maturity, above ground biomass was collected from 1 m of row and dried in forced air dryers at 55 °C for at least 7 days. Total weight, spike count, and spike weight were recorded for all samples. Grain was harvested from manually detached spikes using a laboratory thresher (LD-180, Wintersteiger Inc., Salt Lake City, UT). In 2013, the procedure was modified to measure spike number and grain yield in a weighed subsample (approximately 200 g) of the harvested material. Straw yield per area (SY) was calculated from the difference between total weight and grain yield. Above-ground non-grain biomass (straw + leaves, SNC) was ground on a Wiley 3 mill with a 2 mm screen. N concentration in ground tissue samples was measured by NIR reflectance using a Perten DA7250 spectrometer (Perten Instruments North America, Springfield, IL) calibrated to combustion analysis (LECO FP528, St. Joseph, MI). Grain N concentration (GNC) was determined by NIR reflectance (Perten DA7250 spectrometer) calibrated to combustion analysis (LECO FP528).

Total nitrogen in plants at maturity per area (NTA) was calculated as the sum of grain nitrogen yield (GNY = GNC × GY) plus straw nitrogen yield (SNY = SNC × SY). Nitrogen utilization efficiency (NUtE) was defined as the ratio of GY to NTA. Nitrogen Harvest Index (NHI) was defined as the ratio of N harvested as grain (GY × GNC) to NTA. Grain protein deviation (GPD) was the residual from the linear regression of GPC on GY in each trial. Within each trial, the total N available to the crop was determined empirically as the 95th percentile of NTA (NTAmax) as per Cormier et al. (2013). Nitrogen uptake efficiency (NUpE) at maturity was defined as the ratio NTA/NTAmax. Nitrogen remobilization (NREM) was calculated as the difference between anthesis N content and SNY. Post-anthesis nitrogen uptake (PANU) was calculated as the difference between NTA and anthesis nitrogen content.

Statistical analysis

Data were analyzed by mixed effects analysis of variance using ASReml-R (Butler et al. 2007). Year (Y), and N rate (N), and check cultivar (C) were treated as fixed effects. Replication (β), incomplete block (I), and genotype (g) were treated as random effects, in an extension of models developed by Wolfinger et al. (1997) and Federer (2005).

$$y = \mu + {\text{Y}} + {\text{N}} + ({\text{Y}} \times {\text{N}}) + {\text{C}} + ({\text{Y}} \times {\text{C}}) + ({\text{N}} \times {\text{C}}) + ({\text{Y}} \times {\text{N}} \times {\text{C}}) + \beta :{\text{Y}} + \beta :({\text{Y}} \times {\text{N}}) + \beta :({\text{Y}} \times {\text{N}} \times {\text{C}}) + I + g + (g \times {\text{Y}}) + (g \times {\text{N}}) + (g \times {\text{Y}} \times {\text{N}})$$

Significance of N rate was evaluated using Wald tests. Significance of g, g × Y, g × N, and g × N × Y were evaluated using likelihood ratio tests. Genetic variance (σ 2 g ), genotype × Year variance (\(\sigma_{{g\,\times \,Y}}^{2}\)), genotype × N rate variance (\(\sigma_{{g\,\times \,N}}^{2}\)), genotype × Year × N rate variance (\(\sigma_{{g\,\times \,Y\,\times \,N}}^{2}\)) were summed to estimate total variance associated with genotype, (\(\sigma_{{Total(g)}}^{2}\)), so that the proportion of total genotype-associated variance due to g, g × Y, g × N, and g × N × Y could be calculated for each trait. Broad-sense heritability was calculated on an entry-mean basis as:

$$H^{2} = \frac{{\sigma_{g}^{2} }}{{\sigma_{g}^{2} + \frac{{\sigma_{{g\,x\,Y}}^{2} }}{2} + \,\frac{{\sigma_{{g\,x\,N}}^{2} }}{2} + \,\frac{{\sigma_{{g\,x\,N\,x\,Y}}^{2} }}{4} + \frac{{\sigma_{{error}}^{2} }}{8}}}$$

Because the genotype × year interaction generally contributed a substantial fraction of the variation associated with genotype, for the purposes of the association analyses and the analyses of time-related trends in the germplasm, best linear unbiased predictors (BLUPs) for test genotypes within each year were obtained in a two-step process. First, incomplete-block adjusted estimates for each genotype within each main plot were calculated. Then BLUPs were calculated from a standard split-plot analysis of variance with N rates as main plots and test genotypes as subplots. These genotype BLUPs will be publicly available on the T3 database (http://triticeaetoolbox.org.) following publication of this manuscript. Spearman rank correlations of the genotypes for each trait between the two years of the trial were calculated using the rcorr function in the Hmisc package in R (R Core Team 2016). Genetic improvements in NUE traits were evaluated with simple linear regression models of the NUE trait on release year for the 183 genotypes released after 1960. Based on the significant trends for decreased height and FLO with year of release, and the correlations of height and FLO with the most NUE traits, optimal multiple regression models using the predictor variables of release year, FLO, and height were selected using the Akaike’s Information Criteria (AIC) in the step() function in R. The rate of change with breeding (% yr−1) was calculated as the ratio of the coefficient of release year (slope) to the BLUP of the cultivar, ‘Scout 66,’ which is a long-term check cultivar used in regional cooperative yield testing nurseries.

Genotypic correlations (r g ), phenotypic correlations (r p ) of all traits and the standard deviations of these correlations were estimated in a pair-wise manner from the incomplete-block adjusted estimates for each test genotype in each main plot in each year using multivariate restricted maximum likelihood estimation in PROC MIXED in SAS v9.3 by extension of the method described by Holland (2006). Year, N rate, and replication were modeled as fixed effects; genotype, genotype × year, genotype × N rate, and genotype × year × N rate were modeled as random effects.

Genome-wide association scans (GWAS)

High-density single-nucleotide polymorphism (SNP) marker data from the wheat 92 K iSelect assay (Wang et al. 2014) were used to identify potential marker-trait associations by GWAS. Data for 21,555 SNP were filtered to include only SNPs with a minor allele frequency (MAF) >0.05 and a fraction of lines with missing calls <0.10 to provide a total of 16,052 SNPs. Missing data were imputed with the EM algorithm in rrBLUP (Endelman 2011). Map positions for 14,829 SNP markers were as given in the 90 K consensus map downloaded from http://wheatgenomics.plantpath.ksu.edu on May 22, 2015; 1223 SNP markers were unmapped. SNP assay names are used in the text and tables. SNP data are available online (https://triticeaetoolbox.org/wheat.). Association analysis was conducted with compressed mixed linear model (Zhang et al. 2010) implemented in the GAPIT R package (Lipka et al. 2012) with the option for optimal model selection using the Bayesian Information Criteria. Covariates provided to the model included FLO, HT, and four principal components calculated within GAPIT from the coancestry matrix. The model selection was provided with a kinship matrix calculated as the Realized Relationship Matrix of Endelman and Jannink (2012). A conservative threshold for marker significance was established based on an experiment-wise error rate of 0.05 using the estimated number of effective markers calculated as described by Li and Ji (2005). Linkage disequilibrium (LD) between significant markers was calculated using the LD function in the genetics package (Warnes et al. 2013) in R. Significant marker-trait associations were collapsed to QTLs based on LD relationship: markers with r < 0.8 to the nearest significant marker were treated as distinct QTLs; markers with r ≥ 0.8 to the nearest significant marker were collapsed in QTLs represented by the most highly significant SNP.

Results

Effect of growing conditions

The 2012 and 2013 growing seasons in NE were markedly different: the 2012 season was characterized by an early, warm spring and the 2013 season by a late, cool spring (Guttieri et al. 2015). Crop development consequently was markedly different in the two seasons. In 2012, the Jagger and Settler CL check plots reached anthesis on May-1 and May-6, respectively, and physiological maturity on Jun-5 and Jun-11, respectively. In contrast, in 2013, the Jagger and Settler CL check plots reached anthesis on May-27 and May-31 and physiological maturity on Jun-27 and Jul-1, respectively. Mean grain yields in the 2012 and 2013 trials were 3580 and 5240 kg ha−1, respectively. Mean grain yield in 2012 was below the target yield used to establish N rates (4000 kg ha−1), and the low and moderate N treatments corresponded to 58 and 88% of target N for the mean yield of the trial. Mean grain yield in 2013 was above the target yield used to establish N rates (4700 kg ha−1), and the low and moderate N treatments corresponded to 43 and 63% of target N for the mean yield of the trial.

Effect of nitrogen rate

The additional N provided by the moderate N rate treatment significantly increased grain N yield (9.1%) and grain yield (6.6%, Table 2). The moderate N treatment also significantly increased grain protein deviation, N utilization efficiency and post-anthesis nitrogen uptake. The effects of N rate were similar in the two years of the trial (non-significant N rate × year interaction, data not shown), except for PANU (N rate × year F = 8.1, p < 0.01). Nitrogen rate did not significantly affect flowering date, height, nitrogen harvest index, nitrogen uptake efficiency, or nitrogen remobilization. Genotypes responded similarly to nitrogen rates (non-significant genotype × N rate interaction). The PANU response was complex, with significant genotype × N rate × year interaction.

Table 2 Effect of nitrogen rate, estimated broad-sense heritability, and contributions of genotype (G), genotype × year (G × Y), genotype × N rate (G × N) and genotype × year × N rate (G × Y × N) to the genetic components of variance of flowering date, height, and nitrogen use efficiency traits
Full size table

Genetic variance for nitrogen use efficiency traits

Wheat genotype was a significant source of variation for all NUE traits measured (Table 2) and for flowering date and plant height, and the genotype × year interaction was a significant source of variation for all NUE traits except NREM and PANU (however the three-way interaction of genotype × N rate × year was significant for PANU). The contribution of genotype to the genetic components of variance (Table 2) was less than the contribution of the genotype × year interaction for GY, GNY, NHI, NUpE, and NUtE, which resulted in moderate heritability estimates for these traits (0.25–0.53). Grain protein deviation was the most heritable NUE trait (0.63).

Correlations between traits

The estimates of phenotypic and genotypic correlations of traits from multivariate multi-location analysis (Table 3) demonstrate that plant height and flowering date were positively correlated (r g  = 0.67 ± 0.04; r p  = 0.52 ± 0.03; Table 3). Genotypic correlations of HT and FLO were negative with all NUE traits except NREM, with which they were positively correlated. Earlier and shorter wheat genotypes tended to be more nitrogen use efficient genotypes. NREM and PANU were strongly negatively correlated (r g  = −0.94 ± 0.04; r p  = −0.89 ± 0.01). Genotypic correlations of GNY and GY were negatively with NREM and positively with PANU; the relationships were stronger for GNY than for GY. PANU had relatively weak genotypic correlation (r g  = 0.37 ± 0.12) with NUtE, but strong genotypic correlation (r g  = 0.72 ± 0.10) with NHI. For some combinations of traits, the genotypic and phenotypic correlations were similar (e.g. GY and GNY, GY and NUpE, NUtE and NHI); for many combinations of traits the genotypic correlations had much greater magnitude than the corresponding phenotypic correlations (e.g. NREM and FLO, HT and NUpE, PANU and NHI).

Table 3 Genotypic (above diagonal) and phenotypic correlation (below diagonal) of flowering date (FLO), height, and nitrogen use efficiency traits for n = 299 genotypes analyzed over 2 years and two N treatments
Full size table

Phenotypic data by year

Given the substantial genotype × year interactions, the data for the two years were analyzed separately to generate estimates (BLUPs) of genotype performance within each year for use in subsequent analyses. These data, summarized in Figs. 1 and 2, demonstrate that the difference in early spring temperatures in the two years (2012 was a warm, early spring; 2013 was a cool, late spring) resulted in a 23 day delay in flowering in 2013 relative to 2012. Plants in 2013 were, on average, 18.3 cm taller than in 2012 (Fig. 1). Because of the extended period of pre-anthesis growth in 2013, mean dry matter at anthesis was 11,960 kg ha−1, nearly double the mean dry matter at anthesis (6020 kg ha−1) in 2012 (Fig. 2). Straw yield (SY) at maturity in the two years was not proportional to dry matter at anthesis: SY averaged 6200 kg ha−1 in 2012, but only 9350 kg ha−1 in 2013. Grain N yield in 2013 was 36% greater than in 2012. N remobilization was 143% greater in 2013 than in 2012, while PANU was, on average, positive in 2012 and negative in 2013. On average, less N was measured in above-ground tissues at maturity than at anthesis in 2013. NHI and NUtE were lower, NUpE was higher, and NREM was much higher in the high-yield conditions of 2013 than in 2012. Plants had greater N uptake and remobilization in high-yield conditions, but used the N less efficiently to produce yield and protein.

Fig. 1
figure 1

Distributions of best linear unbiased predictors (BLUPs) of plant height and flowering date for 299 genotypes grown in 2012 and 2013. The center, top, and bottom of the boxes indicate the mean, mean plus one standard deviation, and mean minus one standard deviation, respectively. Red dots indicate the median. Width of the violins are drawn proportional to the frequency distribution of BLUPs. Spearman rank correlation coefficients for genotype BLUPs in the two trial years are indicated for each trait; n = 299 and three asterisks indicates significance of correlation coefficients at p < 0.001

Full size image
Fig. 2
figure 2

Distributions of best linear unbiased predictors (BLUPs) of nitrogen use efficiency traits for 299 genotypes grown in 2012 and 2013. The center, top, and bottom of the boxes indicate the mean, mean plus one standard deviation, and mean minus one standard deviation, respectively. Red dots indicate the median. Width of the violins are drawn proportional to the frequency distribution of BLUPs. Spearman rank correlation coefficients for genotype BLUPs in the two trial years are indicated for each trait; n = 299, ns = non-significant, and three asterisks indicates significance of correlation coefficients at p < 0.001

Full size image

The Spearman rank correlations of the genotypes between the two years of the trial were strongly positive for flowering date (r = 0.68, n = 299, p < 0.001) and plant height (r = 0.79, n = 299, p < 0.001, Fig. 1). Correlations were weaker for grain protein deviation (r = 0.42, n = 299, p < 0.001) and NUtE and PANU (r = 0.35, n = 299, p < 0.001 and r = 0.34, n = 299, p < 0.001, respectively, Fig. 2). Correlations were significant for all traits other than NUpE, and are generally consistent with the rankings of the estimated heritabilities (Fig. 2).

Trends in germplasm

In both 2012 and 2013, simple linear regression of the genotype BLUPs for HT, FLO, and the NUE traits on release year for the 183 cultivars released after 1960 was significant for most traits (Table 4). Anthesis date and NREM regression on release year were weakly significant (p = 0.03 and p = 0.04) in 2012 and non-significant (p > 0.05) in 2013. GPD regression on release year was non-significant (p > 0.05) in 2012, but was highly significant (p = 0.003) in the high-yield conditions of 2013. With the exception of PANU, the slope of the regressions was greater under the high-yield conditions of 2013. The trends in the germplasm indicate increasing PANU, NUpE, NHI, NUtE, GNY, and GY with year of release.

Table 4 Trends in flowering date, height, and nitrogen use efficiency trait trends with year of release for 183 hard winter wheat genotypes released as cultivars since 1960
Full size table

Accounting for FLO and/or HT in multiple regression models with release year substantially improved the accuracy of models for NUE traits as demonstrated by the greater R 2 values (Table 4). In 2012, the trend was for earlier FLO with year of release in the simple linear model. But when HT was accounted for in the regression model, the model fit improved greatly and the trend was for later FLO with year of release. The model coefficients for year of release were attenuated in the multiple regression models, relative to the simple linear models, providing more conservative estimates of genetic progress for NUE traits, accounting for trends in HT and FLO.

Genome-wide association scans (GWAS)

Our objective in GWAS was to identify genomic regions associated specifically with NUE traits, rather than to identify genomic regions associated with height and flowering time. The strong genetic correlations of HT and FLO with NUE traits indicated that HT and FLO would be important covariates in the GWAS. Therefore, a model selection approach with HT, FLO, and principal components supplied as possible covariates was used to identify the optimal model, based on Bayesian Information Criteria, then to evaluate the effects of markers after accounting for these fixed effects. Model selection retained height as a covariate for every NUE trait in each year, and retained FLO as a covariate in all cases, excepting NUpE in 2012. The first principal component was retained as a covariate only for NUtE and GNY in 2013. The covariates, together with the kinship matrix, explained as much as 65% of the variance for individual NUE traits (Table 5).

Table 5 Covariates retained in genome-wide association analyses and peak mapped SNP associated with nitrogen use efficiency traits
Full size table

Applying the liberal threshold of marker significance of p < 0.001, as in Cormier (2015), a total of 224 significant marker associations were identified for the eight NUE traits (GNY, GY, GPD, NHI, NUpE, NUtE, NREM, PANU, Online Resource 2). These associations involved 183 unique SNP markers distributed across every chromosome except 6D and 7D. Sixteen of the SNP markers with significant associations were unmapped. The maximum improvement in model R 2 with the addition of a SNP marker was 0.051. A more conservative threshold for marker significance based on an experiment-wise error rate of 5% and estimating the effective number of marker tests according to Li and Ji (2005) applied a threshold of p < 1.168 × 10−4 for significance. Li and Ji’s method accounts for the correlation among the markers, and estimated that the 16,052 markers effectively provided 439 marker tests. At this threshold, 42 significant marker associations with the eight NUE traits were identified in the two years, 15 in 2012 and 27 in 2013 (Online Resource 2). The marker associations involved 25 SNP markers on 5 chromosomes (1A, 1D, 2B, 2D, and 4B) and three unmapped SNPs. This conservative set of associations will be the focus of the subsequent discussion.

In 2012, IWA693 (Chr 1A, 152.2 cM) was significantly associated with GNY (Table 5). IWA693 also had p < 0.001 association with GY in 2012. Six associations were detected for GPD in 2012, on Chr 1D (112.7 cM), 2B (119.1 cM), 2D (80.1 cM), and three unmapped SNPs. The peak SNP on 2D associated with GPD in 2012 (IWB43945) also was significantly associated with NHI and NUtE in 2012.

In 2013, a SNP on Chr 2D, (IWB74084, 76.6 cM), was significantly associated with GNY, and NUpE (Table 5). This SNP was in weak LD (r = 0.30) with IWB43945, which was important in 2012. Two SNPs on Chr 4B, (IWB27326 and IWA7566, located at 61.8 and 80.6 cM, respectively), were significantly associated with GPD in 2013. These SNPs were in higher LD (r = 0.50) than would be expected based on the consensus map positions. Both SNPs had low minor allele frequency (MAF), ≤ 17%. A set of 10 SNPs on Chr 4B that were in high linkage disequilibrium with each other (r > 0.9), represented by IWB34983 (64.26 cM), were significantly associated with NHI in 2013. Two additional SNPs on Chr 4B in moderately high LD with this block of SNPs (r > 0.7) also were significantly associated with NHI. IWB34983 and SNPs in high LD with IWB34983 also were significantly associated with NUtE in 2013 and had p < 0.001 association with GNY and GY in 2013. No significant marker associations with GY, NREM, or PANU were detected in either 2012 or 2013 at the more conservative threshold of p < 1.168 × 10−4.

The minor alleles of the most significant SNPs generally had adverse effects on NUE traits (Table 5), which indicates that favorable alleles for NUE have been selected in breeding. Exceptions, cases in which the less common allele was favorable, were the 2B and 2D SNPs for GPD in 2012, and the 2D SNP associated with GNY and NUpE in 2013. Increasing the frequency of these alleles in hard winter wheat germplasm may lead to improved NUE.

Discussion

Genetic variance for NUE traits

Genetic variation for NUE traits is present within the Great Plains hard winter wheat germplasm pool. Among the eight NUE traits evaluated in this study, grain protein deviation, post-anthesis nitrogen uptake, and N utilization efficiency had the greatest broad sense heritability. Genotype × environment interaction effects were highly significant, likely due to variation in genetic responses to the very different spring temperature conditions in the two years of the trial. The proportion of genetic variance for GY due to genotype (35%) and genotype × year interaction (62%) were substantially different than in the European study reported by Cormier et al. (2013) (genotype = 60%, genotype × env = 36%). Mean grain yield (4670 kg ha−1) and the heritability of grain yield (0.43) were substantially lower in our study than in European trials (mean grain yield = 7400 kg ha−1, heritability = 0.79). However, GNY was nearly identical to the European study, due to low grain protein concentration (mean 9.93%) in Cormier et al. (2013), compared to 15.2% and 15.9% in our study in 2012 and 2013, respectively. Nitrogen uptake efficiency was similar to the European study, but nitrogen utilization efficiency (22.9) was lower than in the European study (48.8). Nitrogen harvest index was greater (0.812) in the European study than in our work (0.615). Genotype × N interaction was not an important source of variation in our study, perhaps because the two rates of N were similar. The high grain mean protein concentrations in our study (>15%) indicate that the crop was not grown under N-limiting conditions.

Reducing post-anthesis N loss, and improving PANU, would be tremendously valuable to NUE. We observed significant genetic variation for the trait, although there was a complex relationship of genotype, N rate, and year. Mean PANU among the genotypes in 2013 was negative (Fig. 2, −66.4 kg N ha−1), while in 2012, mean PANU was positive (19.7 kg N ha−1). Post-anthesis nitrogen loss of similar magnitude has been reported previously for wheat (Papakosta and Gagianas 1991). Ammonia loss from the tops of plants is reportedly the major mechanism of N loss (Harper et al. 1987; Farquhar et al. 1980), and losses are greatest under conditions of high evapotranspiration rates: low humidity, high temperature, and good available soil moisture. In the 14 days following mean anthesis date, the average daily temperature was 18.8 C in 2012 and 18.7 C in 2013. Average relative humidity was somewhat lower in this period during 2013 (73.6%) than in 2012 (81.7%). Soil moisture conditions following anthesis were very good in both years. Therefore, environmental conditions following anthesis may not explain the difference between the two years in PANU. The mean anthesis date in 2012 was 23 days earlier than in 2013 (Fig. 1), which shortened the growing season and contributed to reduced average total dry matter at anthesis in 2012. This difference in dry matter at anthesis in the two years was reflected in lower N content at anthesis, 144 kg N ha−1, in 2012, compared to 317 kg N ha−1 in 2013. Yet average grain N yield was only 36% greater in 2013 than in 2012, and grain yield was only 29% greater in 2013 than in 2012 (Fig. 2), which suggests that sink capacity may not have increased in proportion to source in the long season, high yield conditions in 2013. After accounting for the effects of plant height and flowering time, no significant SNP associations with PANU were identified at the conservative threshold, and only 11 significant SNP associations on Chr 5A, 2B, 2D, 3D, and 5D were identified at the more liberal p < 0.001 threshold (Online Resource 2).

Genotypic correlations among traits

Correlations among NUE traits are expected based on the numerical interdependence of the traits. Genotypic and phenotypic correlations of traits (Table 3) were most similar for traits with numerical interdependence. For example, the nitrogen utilization efficiency trait uses grain yield as the numerator in the calculation, and nitrogen harvest index uses grain nitrogen yield as the numerator in calculation. PANU and NREM were strongly negatively correlated (r g  = −0.94 ± 0.04; r p  = −0.89 ± 0.01). This correlation is based on the calculation of the traits. PANU is the difference between total N yield at maturity and total N yield at anthesis. As N yield at anthesis increases, other terms being equal, PANU decreases and NREM increases; and as straw N at maturity increases, other terms being equal, PANU increases and NREM decreases. Among the NUE traits, GPD generally was least correlated with the other traits. This would be anticipated because it is a measure of deviation from the physiological relationship between yield and grain protein. The absence of strong genetic correlations with other traits suggests that GPD may be manipulated independently.

In contrast, genetic correlations of plant height and flowering date, and genetic correlations of these traits with NUE traits, indicate pleiotropy. Flowering date and plant height genetic correlations with NUE traits generally were very strong, while phenotypic correlations with NUE traits generally were much weaker (Table 3). Phenotypic correlation of two traits, X and Y, results from the combined effect of both genetic and environmental correlation:

$$r_{p} = h_{X} h_{Y} r_{g} + \,e_{X} e_{Y} r_{E}$$

where h is the square root of the heritability, e is the square root of the proportion of phenotypic variance attributable to environment, E (Falconer and Mackay 1996). When heritability is low, and when r E is low, phenotypic correlation may be much less than genotypic correlation. In a review of 32 plant studies in which both phenotypic and genotypic correlations were reported, Waitt and Levin (1998) found that the magnitude of correlation in genotypic correlation matrices was greater than in phenotypic correlation matrices in 85% of studies. They also noted that the traits for which phenotypic and genotypic correlations were most similar for traits that belonged to the same functionally and/or developmentally related set that they described as being “a higher degree of genetic and phenotypic character integration.” The important effect of flowering time has been observed in phenotypic correlations reported in other NUE studies (Bogard et al. 2011a; Laperche et al. 2007; Cormier et al. 2013). The magnitude of the genotypic correlations of FLO and HT with NUE traits in our study were much greater than the magnitude of phenotypic correlations reported previously.

Genetic improvement in NUE traits

Wheat breeding within the Great Plains of the United States has improved nitrogen use efficiency: within the subset of 183 cultivars released from 1960 to 2014, coefficients of the regression on release year for key NUE traits were significant and positive in both years of the study (Table 4). Accounting for trends in flowering time and plant height attenuated but did not negate the estimates of gain in NUE parameters. Genetic improvement of NUtE was positive in both 2012 (0.115% yr−1) and 2013 (0.367% yr−1), after accounting for changes in HT and FLO (Table 4). However, the rate of genetic improvement we observed for NUtE was approximately one-half the rate of improvement we observed for grain yield in 2012 (0.331% yr−1) and 2013 (0.761% yr−1). Nitrogen harvest index improved at a rate of 0.05% yr−1 in 2012 and 0.20% yr−1 in 2013. Previous estimates for the rate of improvement of NHI in European germplasm were 0.15% yr−1 (Brancourt-Hulmel et al. 2003) and 0.12% yr−1 (Cormier et al. 2013). It appears, therefore, that yield gains have come primarily through increased dry matter, rather than through efficient uptake and remobilization of N.

Association analysis

At a liberal threshold of marker significance, a total of 224 significant marker associations involving 183 unique SNP markers were identified for the eight NUE traits. These markers were distributed across every chromosome except 6D and 7D. QTL effects were modest: the maximum improvement in model R 2 with the addition of a SNP marker was 0.051. At a more conservative threshold for marker significance, a threshold of p < 1.168 × 10−4 for significance, 42 significant marker associations involving 25 SNP markers on 5 chromosomes (1A, 1D, 2B, 2D, and 4B) and three unmapped SNPs were identified in the two years. Previous QTL studies have identified between 32 and 380 QTLs for NUE traits (Cormier et al. 2016).

The strong genotypic correlations of flowering date and height with NUE traits was reflected in the high model R 2 for most of the trait-year combinations in the absence of any fitted SNPs. Despite attempting to control for the effects of plant height in the GWAS models, SNP markers on Chr 4B were significantly associated with GPD and NHI in 2013. Rht-B1b, located on chromosome 4B, became the prevalent allele for semi-dwarf stature in modern Great Plains hard winter wheat after the late 1970s (Guedira et al. 2010). Grogan et al. (2016) found that 72% of the genotypes included in this study were homozygous for the semi-dwarf-allele Rht-B1b. The significant SNP markers on 4B in this study were in moderate LD (r = 0.58) with the Rht-B1 marker data previously reported for 291 of these genotypes (Grogan et al. 2016). The GPD-associated SNP alleles are at a relatively low frequency in the population of lines (14–17%), and have a negative effect on GPD, while the NHI-associated SNP alleles are present at a relatively high frequency (47%) and likewise have a negative effect on NHI and on NUtE. Therefore the favorable 4B alleles for NUE traits appear to be widely distributed in the Great Plains hard winter wheat germplasm. The 4B QTLs may not be a consequence of Rht-B1. GSe, a cytosolic gene encoding glutamine synthetase, maps near Rht-B1, and the temporal profile of its activity suggests that it is a key gene involved in N remobilization (Habash et al. 2007; Fontaine et al. 2009).

In both years, SNP markers on the long arm of Chr 2D were associated with NUE traits, which suggests that this may be a fertile target for future investigation. The favorable 2D alleles for NUE traits were present at low frequency (0.13–0.18) in the germplasm included in this study. Previous studies that have identified NUE associations with 2D have localized on or near the Ppd-D1 locus for photoperiod sensitivity on 2DS (Cormier et al. 2014; Laperche et al. 2007). In contrast, the SNPs we have identified map to 2DL and are not in LD with previously reported Ppd-D1 marker data for these genotypes (Grogan et al. 2016). Glutamine synthase 2 genes are located in distal deletion bins on the long arms of group 2 chromosomes and are potential candidates for future investigation, given their important role in previous studies (Li et al. 2011).

Conclusions

Nitrogen use efficiency has improved in the Great Plains hard winter wheat germplasm in the period from 1960 to 2014, even after controlling for the effects of changes in height and heading date in this time period. Genotypic correlations of height and flowering date with NUE traits other than grain protein deviation were highly significant, and height and flowering date were retained as covariates in GWAS models for NUE traits in both years of the trial for most traits. In many cases, height and flowering date explained substantial proportions of the variation in the GWAS. Including height and flowering date as covariates may be a useful approach in other studies. Our results point to the difficulty inherent in identifying alleles for improved NUE within un-adapted germplasm: confounding effects of non-adaptive phenology may obscure modest improvements in N uptake, assimilation, or remobilization. Moreover, bi-parental populations with uniform phenology may be most effective for mapping of NUE QTLs.

Using a conservative threshold of marker significance for GWAS, a limited number of significant marker-trait associations were identified in the Great Plains hard winter wheat germplasm that may have future utility in breeding. However, relative to the effects of flowering time and height, marker effects were small. Consistent, favorable associations with NUE traits were identified on 2DL that merit further investigation within germplasm with uniform flowering time and plant height to better quantify the marker allele effects independent of phenological adaptation. The favorable alleles on 2DL were relatively rare in the Great Plains hard winter wheat germplasm, indicating that 2DL may be a fertile region for continued study.