Many wheat (Triticum aestivum L.) production regions are threatened annually by drought stress. Carbon isotope discrimination (Δ) has been identified as a potentially useful trait in breeding for improved drought tolerance in certain environments. Broad use of Δ as a selection criterion is limited, however, mainly due to an inconsistent relationship observed between grain yield and Δ and, to a lesser extent, because of the high resource demand associated with phenotyping. The efficiency of selection may be improved by the identification and verification of molecular markers for use in marker-assisted selection (MAS), and a reliable relationship to grain yield may be established based on a location’s total amount and distribution of precipitation over the growing season. Given the environmental variability in precipitation dynamics, it is necessary to evaluate this relationship in target breeding environments. In this study, grain Δ was collected on a panel of 480 advanced soft white winter wheat varieties grown in five Pacific Northwest environments. A genome-wide association study approach was used to evaluate the amenability of grain Δ to MAS. The genetic architecture of grain Δ was determined to be characterized by multiple, small effect marker-trait associations with limited repeatability across environments, suggesting that MAS will be ineffective at improving Δ selection efficiency. Further, the relationship between grain yield and Δ ranged from neutral () to moderately positive () in the target environments. Such moderate correlations, coupled with variability in this relationship, indicate that direct selection for Δ may not be beneficial.
Drought stress is the single greatest threat to wheat (Triticum aestivum L.) productivity worldwide, impacting more than half of all growing regions [1, 2]. As one of Washington State’s most widely cultivated crops, wheat and the associated industries are critical to rural economies and to the State economy as a whole . Like much of the world, drought represents one area of significant vulnerability to Washington’s wheat growers . In 2015, drought stress resulted in an estimated 22 percent reduction in wheat yield, representing a total economic loss of $212.4 million . Though an extreme case, this was not an isolated event. Washington’s major wheat-producing area is among the most drought-prone regions in the country .
Plant breeding efforts to develop more tolerant varieties of wheat may be able to mitigate the impacts of drought. Breeding for drought tolerance is complicated, however, largely due to the quantitative nature and significant genotype-by-environment interactions associated with drought tolerance. Past efforts targeting drought tolerance directly based on yield components have experienced limited success . An alternative strategy, involving selection for less complex physiological traits that are associated with yield in dry environments (leaf rolling, stay green, canopy temperature depression, osmotic adjustment, water soluble carbohydrates, deep roots, and glaucousness), has been suggested as a more effective approach . Many of these traits are less influenced by the environment and can therefore contribute to greater yield stability .
Carbon isotope discrimination (Δ) has been identified as an indirect selection criterion for improving grain yield in certain environments [10, 11]. In the early 1980s, researchers observed that the molar concentration ratio of 13C relative to 12C in plant tissue is often lower than the same ratio of atmospheric carbon. This finding indicates that plants discriminate in favor of 12C during the conversion of CO2 into biomass. As a standard for evaluating plant isotopic composition, was proposed by Farquhar and Richards  and is defined as the deviation from unity of atmospheric 13C/12C over plant 13C/12C.
Wide genetic variation in Δ exists among and between plant species. This variation may be largely explained by genotypic differences in stomatal conductance and photosynthetic capacity. Early research found Δ to share a positive and linear relationship with the ratio of internal leaf CO2 concentration () and the atmospheric CO2 concentration () . This ratio is principally impacted by stomatal conductance and photosynthetic capacity . A high Δ genotype may therefore be assumed to express high stomatal conductance, low photosynthetic capacity, or both, whereas the opposite may be assumed for a low Δ genotype. The relationship to makes it possible to draw inferences regarding fundamental plant physiology based on measurements of Δ. Further, Δ can be highly heritable and relatively easy to manipulate in breeding programs [15, 16].
Despite these positive characteristics, common use of Δ as a selection criterion has been limited. This is largely due to the variable relationship observed between Δ and grain yield across environments and, to a lesser extent, because of the resources necessary to phenotype. Selection efficiency may be improved by the identification and verification of molecular markers associated with Δ for use in marker-assisted selection (MAS). A genome-wide association study (GWAS) can be considered one of the first steps in identifying potential candidate markers for MAS . By way of statistical analysis, a GWAS detects DNA markers that share an association to the phenotype [18, 19]. The results provide researchers with a better understanding of the trait’s genetic architecture, including the number of marker-trait associations (MTAs) involved in trait expression as well as the relative contribution of each genetic region to the observed phenotype. For effective use in MAS, the identified DNA markers should explain a high proportion—typically between 40 and 60%—of the trait’s phenotypic variation and demonstrate repeatability across environments [11, 20, 21]. Although Δ is influenced by both stomatal conductance and photosynthetic capacity, traits that are polygenic themselves, there is evidence for moderate effect markers, explaining 30% of the phenotypic variance , as well as marker repeatability . These findings suggest that Δ may be amenable to MAS, warranting further examination.
Reasons for the inconsistent grain yield and Δ relationship are not well understood, though differences in environmental conditions, principally the distribution and amount of precipitation, may be large factors. Most field studies conducted under moderate water limitations report either a positive or neutral relationship between grain yield and Δ [24–28]. Only infrequently are negative correlations observed, occurring most often in very dry environments where crops experience little precipitation during the growing season and rely on stored soil moisture as a water source [10, 15, 29]. For effective use as a selection criterion, the relationship to grain yield in a particular environment ought to be consistent from year-to-year, either positive or negative, and of reasonably high correlation (typically or ). The complexity between Δ and grain yield reported in the literature, along with the potential to improve yield under drought stress, makes it appropriate to study this relationship in target environments where breeding lines are tested.
The relationship between Δ and grain yield has not been evaluated for the primary wheat-producing region of Washington State. Precipitation across this region follows a low (<300 mm) to high (450-600 mm) gradient moving east from the foothills of the Cascade Mountains to the Washington-Idaho border . The winter wheat breeding program at Washington State University targets two main environments—low and high precipitation. The objective of this study was to explore the relationships between Δ and grain yield and, more broadly, to evaluate the utility of Δ as a selection criterion to target breeding environments in Washington State. To this end, practical elements of selection are considered beyond the trait’s relationship to yield, including trait heritability, genetic variation, and ease of selection.
2. Materials and Methods
2.1. Plant Material
This study collected data on a panel of 480 advanced soft white winter wheat varieties from U.S. Pacific Northwest breeding programs (Oregon State University, University of Idaho, Washington State University, USDA-ARS, and private breeding companies). All varieties included in the panel were from the soft white market class, of which 63% were lax-headed, and the remaining 37% were club head type. From the 480 lines in the original panel, a total of 459 lines were ultimately included in the association analysis. The 21 excluded lines were either identified as a separate phenotypic class, had poor DNA quality, or were missing marker covariate information. Previous work on this panel established the genotype data for each variety . The panel was genotyped using an Illumina Infinium iSelect 90K single nucleotide polymorphism (SNP) chip at the USDA-ARS Biosciences Research Laboratory in Fargo, ND, USA.
2.2. Experimental Design
A panel of 480 advanced soft white winter wheat varieties was grown in an unreplicated augmented block design in five field environments. The panel was phenotyped for plant height (cm), heading date (Julian), grain yield (t ha-1), and Δ. Location-specific monthly precipitation throughout the growing season was obtained from databases of nearby weather stations. A GWAS was conducted on the panel for grain yield and Δ. MTAs determined to be significant were analyzed for physical linkage, repeatability across environments, percentage of phenotypic variation explained (2), and colocalization with significant grain yield MTAs. Summary statistics, as well as inter- and intraenvironment phenotypic variance, for each trait was calculated. Pearson correlations between Δ and grain yield, plant height, and heading data were also determined (Figure 1).
2.3. Field Conditions
The panel was grown in an augmented block design with unreplicated 2.8 m2 plots at the Spillman Agronomy Farm near Pullman, WA (46°7N; -117°1W), during 2015, 2016, and 2017. In 2017, the panel was also grown in an augmented block design with unreplicated 6.1 m2 plots near Lind, WA (46°8N; -118°6W), and 8.7 m2 plots near Pendleton, OR (45°7N; -118°6W). In each location-year, 20% of the plots were planted to the check cultivar “Madsen” (PI 511679) . Trials were planted in Pullman on October 8, 2014, October 8, 2015, and November 3, 2016; in Lind on September 7, 2016; and in Pendleton on October 4, 2016. Trials were harvested in Pullman on July 23, 2015, August 5, 2016, and August 10, 2017; in Lind on July 20, 2017; and in Pendleton on July 25, 2017. The mean annual precipitation for these locations (2010-2017) is 230 mm at Lind, 435 mm at Pullman, and 235 mm at Pendleton. The mean annual temperature (2010-2017) for Pullman is 8.8°C. In Lind and Pendleton, the mean annual temperature (2010-2017) is 10.2°C and 11.8°C, respectively. Precipitation and temperature data for each location-year were collected from weather stations near the field sites. AgWeatherNet sites (https://weather.wsu.edu/) were used to access precipitation data for all Pullman environments, as well as Lind. An AgriMet weather station (https://www.usbr.gov/pn/agrimet/wxdata.html) entitled “Echo” was used to access precipitation data for the Pendleton field site. This station is approximately 38 miles from the Pendleton field location, whereas the two AgWeatherNet stations are at the same location as the Pullman and Lind field sites.
The amount of precipitation that each location-year received during the growing season ranged from low in Pendleton (265 mm) and Lind (294 mm) to high in Pullman 2016 (463 mm). Moderate levels of precipitation were experienced in Pullman 2015 (384 mm) and Pullman 2017 (377 mm) (Table S1). Both Lind and Pendleton had fairly consistent precipitation throughout the vegetative stages. Pullman 2016 and Pullman 2015 experienced a general pattern of declining precipitation after the seedling stage. Pullman 2017 received high early precipitation and high precipitation in March, with relatively little precipitation during December and January. All environments tended to experience a decrease in precipitation from spring into summer. Pullman 2015 received the least postanthesis precipitation (10 mm) followed by Lind (14 mm), Pullman 2017 (21 mm), and Pullman 2016 (32 mm). Although heading date was not collected for Pendleton 2017, the final two months of the growing cycle received 19.3 mm of precipitation.
2.4. Phenotypic Data
Plot heading date (Julian) was determined as the number of days from January 1 to 50% of fully exposed heads. Heading date was recorded in all trials except in Pendleton 2017. Plant height (cm) was measured from the base of the plant to the top of fully emerged heads, excluding the awns. Grain yield (t ha-1) was determined from the grain weight per plot obtained from a ZÜRN 150 plot combine harvester (ZÜRN Harvesting, Germany). Grain samples (20 g) were collected from harvested seed for each variety and milled into flour using an UDY cyclone sample mill (UDY Corporation, CO). Flour samples weighing between 3.00 and 5.00 mg were packed into sterile tin capsules. Grain carbon isotope percentages were measured using an isotopic mass spectrometer at the Stable Isotope Core Laboratory at Washington State University. Carbon isotope composition (δ13C) was calculated by comparing the ratio of 13C to 12C for each sample (s) against the same ratio of a Vienna Pee Dee Belemnite (VPDB) standard (VPDB) by the following formula: . The Δ value of the sample was then calculated as , where δ13Ca and δ13Cp are the carbon isotope compositions of the atmosphere and plant samples, respectively . The δ13Ca value was assumed to be −8.0 per mil . Carbon isotope composition values are unitless, though they are often expressed in terms of “per mil,” indicating that the original value was multiplied by 103. Doing so is a matter of convenience, allowing values to appear as whole numbers.
2.5. Statistical Analysis
The phenotypic data were assessed for influential outliers (greater than three standard deviations). Three suspicious data points were identified in Pullman 2015, six in Pullman 2016, four in Pullman 2017, two in Lind 2017, and four in Pendleton 2017. These outliers were excluded from any further statistical analyses. Pearson correlations and summary statistics (mean, range, and standard deviation) were calculated for phenotypic traits using the R statistical software (R Core Team, 2016). Given the potential to confound the relationship between grain yield and Δ, a separate correlation analysis was conducted between grain yield and Δ excluding those genotypes with heading date greater than one standard deviation and less than negative one standard deviation from the environment-specific mean. Variation between environments and between genotypes and checks within environments was calculated using SAS version 9.4 (SAS Institute Inc., Cary, NC, USA). ANOVA coupled with post hoc Tukey HSD and linear contrast tests calculated in R statistical software were used to identify any significant differences between environments for phenotypic traits.
Broad-sense heritability for each trait was determined based on the formula , where is the genetic variance and is the phenotypic variance . The for each trait was calculated by dividing the sum of squared residuals of an equal means model by the number of observations. The environmental variance, , was calculated as the sum of the squared residuals for each genotype across environments divided by the total observations. The relationship allowed for an estimate of and subsequently an estimate of for each trait.
Best linear unbiased predictions (BLUPs) were calculated for Δ and grain yield using the lme4 package in R. The mixed linear model considered genotype and environment as random effects. Due to the potential to confound Δ results, both heading date and plant height were added to the Δ mixed linear model as covariates. The proportion of phenotypic variance explained by significant SNPs () identified in each location-year and BLUP analysis was calculated in JMP v. 8.1 by fitting genotype and phenotype data on a stepwise regression model and computing the difference in variance explained between a model including all significant SNPs and a model excluding the marker under consideration.
2.6. Linkage Disequilibrium (LD)
Linked SNPs determined to be significantly associated with the phenotype were evaluated for physical linkage. The intersection of the locally weighted polynomial regression- (LOESS-) based curve of the LD parameter and the 95th percentile of the unlinked distribution was considered to be the critical threshold beyond which LD was due to physical linkage . Previous research on the panel used for this study determined that SNPs with an value greater than 0.18 and within 4 cM could be considered to be on the same linkage block . The pairwise LD parameter was computed in TASSEL version 5.2.25 .
2.7. Association Analysis
Association analysis for grain yield and Δ was performed on both the location-year phenotypic data and the BLUP data using 15,229 high-quality SNP markers. Those markers with more than 20% missing data or with a minor allele frequency (MAF) less than 0.05 were excluded from the association analysis. An assessment of the panel’s population structure was previously conducted . In the grain yield and Δ association analyses, model fit was verified by evaluating quantile-quantile plots (Figures S1 and S2). In the Δ association analysis, principal components 1 and 2 were included in the final model to account for the population structure. If heading date or plant height was significantly () correlated with Δ in a particular environment, the covariate was added to the final model for that environment. Heading date was added as a covariate to Pullman 2015 and Lind 2017. Plant height was added as a covariate to all environments excluding Lind 2017. The BLUP analysis included both plant height and heading date as covariates. For the grain yield association analysis, model fit did not improve with the addition of principal components, so a model with no covariates was ultimately selected. The purpose of this analysis was to identify significant grain yield MTAs that colocalize with significant Δ MTAs, making a basic model for grain yield association analysis appropriate. That said, in addition to a colocalization analysis with no covariates in the grain yield model, the same analysis was conducted using a grain yield model including principal components 1 and 2. No notable differences in the regions or number of MTAs found to colocalize were identified between the two analyses.
Fixed and random Circulating Probability Unification (FarmCPU), implemented in GAPIT2 , was used for the marker-trait association analysis. Associations were determined to be significant by the Bonferroni correction method  with —equivalent to a marker-wise threshold value of (0.05/15,229 SNPs). As the Bonferroni correction is overly conservative, assuming that each association test is independent of all other tests , the colocalization assessment considered SNPs significant with a marker-wise threshold value of .
3.1. Phenotypic Variance of Measured Traits and Heritability
A wide variation was observed for all measured traits in each environment (Table 1). The mean grain yield varied from 3.0 t ha-1 in the Lind 2017 trial to 7.8 t ha-1 in the Pullman 2017 trial. The variation for grain yield was significant between the five environments (), as well as between genotypes within each environment (). A post hoc analysis revealed that all five environments were significantly different for grain yield. Mean Δ was lowest in the Lind 2017 trial (16.5) and highest in the Pendleton 2017 trial (18.2). The range of Δ values within environments was roughly 2.0 per mil for all environments excluding Pullman 2015, which had a range of 3.1 per mil. The variation in mean Δ values between environments was significant (), as was the variation in Δ values between genotypes within each environment (). A post hoc analysis revealed that Pullman 2017 and Pendleton 2017 were the only environments without a significant difference for Δ values (). In addition to grain yield and Δ, a significant variation for plant height and heading date was observed both between environments () and between genotypes within each environment (). Considering the Madsen check cultivar only, there were no significant differences for Δ (), grain yield (), plant height (), or heading date () within each environment, providing evidence for limited spatial field variation. Based on this evidence, no data transformations or adjustments were conducted. Broad-sense heritability for Δ considering all environments together was determined to be 0.15, whereas broad-sense heritability for grain yield was 0.06 (Table S2). Broad-sense heritability of Δ calculated for the Pullman location-years separately increased to 0.27.
3.2. Correlation Analysis of Agronomic Traits and Δ
Grain yield and Δ were positively correlated () in all environments, excluding Pullman 2015, where no correlation was observed (Table 2; Figure 2). For those environments with a positive relationship, Pearson phenotypic correlations ranged from 0.09 in Lind 2017 to 0.44 in Pullman 2016. Plant height and heading date were both significantly () correlated with Δ in all environments for which the data was collected. The Pearson correlations between plant height and Δ ranged from -0.13 in the Lind 2017 trial to -0.49 in the Pullman 2017 trial. The range of Pearson phenotypic correlations between heading date and Δ varied from -0.09 in Pullman 2016 to -0.25 in Pullman 2015. In a separate analysis, in which early and late heading genotypes were excluded, the correlation between grain yield and Δ for each environment remained similar to the original correlation analysis (Table 2).
3.3. Model Selection for Association Analysis
Model fit for the grain yield and Δ association analyses was verified with quantile-quantile plots (Figures S1 and S2). In the Δ association analysis, principal components 1 and 2 were included in the final model to account for the population structure. If heading date or plant height was significantly () correlated with Δ in a particular environment, the covariate was added to the final model for that environment. Heading date was added as a covariate to Pullman 2015 and Lind 2017. Plant height was added as a covariate to all environments excluding Lind 2017. The BLUP analysis included both plant height and heading date as covariates. For the grain yield association analysis, model fit did not improve with the addition of principal components, so a model with no covariates was ultimately selected. The purpose of this analysis was to identify significant grain yield MTAs that colocalize with significant Δ MTAs, making a basic model for grain yield association analysis appropriate. That said, in addition to a colocalization analysis with no covariates in the grain yield model, the same analysis was conducted using a grain yield model including principal components 1 and 2. No notable differences in the regions or number of MTAs found to colocalize were identified between the two analyses.
3.4. Analysis of Significant MTAs
A total of 25 MTAs were significant () for grain Δ across 13 chromosomes (Table 3). The number of significant MTAs varied by environment ranging from two identified in Pendleton 2017 to six in Lind 2017. There were five significant MTAs in the Pullman 2015 environment. In Pullman 2016 and Pendleton 2017 and in the BLUP analysis, four MTAs were found to be significant. The majority of these MTAs were detected on the B genome (14), followed by the A genome (8) and the D genome (3). The most informative chromosomes were 4B, with four MTAs, and 5B, with six MTAs. None of the identified MTAs explained greater than 10% of the observed phenotypic variation (), and only four explained greater than 5%.
Four MTAs and three linkage blocks were found to be significant () in more than one environment (Tables S3 and S4). However, none of these regions were significant in more than two environments. Two repeated regions were identified on chromosome 4B. On chromosome 5B, two repeated MTAs and one repeated region were identified. Significant () MTAs for Δ identified in a particular environment were assessed for colocalization with significant () MTAs for grain yield identified in that same environment (Table S5). One linkage block on chromosome 6B was found to be significant () for both Δ and grain yield in the Lind 2017 environment. SNP IWB56415 (60.8 cM) was significant for Δ in Lind 2017, and SNP IWB29847 was significant for grain yield (64.7 cM). The pairwise LD value between these SNPs was 0.24. In the BLUP analysis, a second linkage block common to Δ and grain yield was identified on chromosome 7B (). SNP IWB5954 (171.11 cM) was found to associate with Δ, and SNP IWB12006 (171.11 cM) was significant for grain yield. The pairwise LD value between these SNPs was 0.85.
4.1. Correlation Analysis of Agronomic Traits and Δ
This study found Δ to be positively correlated with grain yield in four of the tested environments and neutral in one of the tested environments (Figure 2). These results are consistent with positive correlations reported in previous studies in which grain Δ was measured from wheat grown under moderate drought conditions [24–28]. Various hypotheses explain this association between grain yield and Δ. High Δ genotypes may achieve high productivity as a result of deep roots that allow access to water . Alternatively, stem reserve carbohydrates, fixed under well-watered conditions, can elevate the overall grain Δ signature once translocated from stem tissue to grain tissue . High Δ grain tissue can also result from stomatal insensitivity to drought, which may be a beneficial trait to a moderate terminal drought environment . Finally, high yield from high Δ genotypes may be the result of an early heading date, allowing the genotype to avoid terminal drought conditions.
Variation in heading date may have been partly responsible for the positive correlations observed in this study. Heading date was significantly associated with Δ in every tested environment. The broadest range in heading date among genotypes spanned 18 days in Pullman 2017, whereas the narrowest range was 11 days in Pullman 2015. The water availability between early and late heading genotypes was likely significant and undoubtedly impacted grain Δ measurements. That said, after removing those genotypes with heading date greater than one standard deviation and less than negative one standard deviation from the mean in each environment, very little difference in the Pearson correlation was observed between grain yield and Δ (Table 2), suggesting that the impact of heading date on the relationship between Δ and yield may be a minor factor in this study.
A number of studies report a repeated and high correlation value between grain yield and Δ, leading researchers to propose grain Δ as a predictive selection criterion for wheat grain yield [10, 42–45]. In the present study, such use of grain Δ cannot be recommended, as the observed correlation was generally moderate. Further, the relationship between grain yield and Δ was inconsistent within the two target breeding environments—high (450-600 mm) and low (<300 mm) precipitation. The high precipitation environment of Pullman saw a neutral relationship in 2015 (; 384 mm precipitation) and a moderately positive relationship in 2017 (; 377 mm precipitation), and in Pullman 2016, the relationship was the most positive of all environments (; 463 mm precipitation). As for the low precipitation environments, Lind 2017 had a slightly positive correlation (; 294 mm precipitation) while the relationship increased in Pendleton 2017 (; 265 mm precipitation). This variation within target breeding environments, combined with the generally moderate correlation observed between grain yield and Δ, makes it difficult to recommend direct selection of grain Δ in either high or low precipitation environments.
4.2. Analysis of Significant MTAs
The genetic architecture of Δ is complex. This appears to be so for Δ measured on unstressed plant tissue, as well as tissue exposed to drought stress. In the case of unstressed tissue, the Δ signature may reflect genetic factors of both stomatal conductance and photosynthetic capacity. These traits are themselves polygenic, and so, a complex genetic architecture may be expected when combined in Δ analysis . Fittingly, a total of 54 QTL have been reported for Δ across a number of mapping studies . In a biparental mapping study of Δ from unstressed leaf tissue, Rebetzke et al.  report Δ QTL on chromosomes 1A, 1D, 2A, 2B, 2D, 3B, 4A, 4B, 4D, 5A, 5B, 6B, 6D, 7A, and 7B. The identified QTL were of small genetic effect, each explaining less than 10% of the additive genetic variance . In the only other known association analysis where Δ is considered, Mora et al.  report significant associations of grain Δ on chromosomes 1A, 3A, 4A, 4B, 4D, 5A, 5B, 5D, 6B, and 7D, with maximum explained phenotypic variation of 9.7% .
In the event that grain fill coincides with a moderate to severe drought, the genetic architecture of grain Δ is likely to be even more complex. Along with genetic contributions of photosynthetic capacity and stomatal conductance, grain Δ may also be impacted by genetic factors indirectly associated with stomatal conductance, including flowering time, root characteristics, and reliance on stem reserve carbohydrates for grain fill . The grain Δ signature depends mainly on the stomatal conductance and photosynthetic capacity occurring at the time of carbon assimilation . Exposure to drought tends to reduce stomatal conductance, driving Δ values down . Those genotypes able to maintain open stomata by accessing water deep in the soil profile likely register higher Δ values. Similarly, stem reserves laid down prior to the drought onset may also register as high Δ, provided fixation occurred during a period of relatively high stomatal conductance. Upon translocation to the grain, these high Δ carbohydrates will elevate the Δ signature of the kernel tissue . Alternatively, high grain Δ could result from early flowering . As plants mature into the later months of summer, it is typical of Mediterranean climates for precipitation and humidity to decline and for temperature to increase. These changes elevate the threat of drought and the likelihood of stomatal closure. A genotype with an earlier heading date may simply be high for Δ by avoiding a period of heightened drought potential.
The relationship that Δ has with makes it possible to approximate the severity of terminal drought in a particular environment. In a well-watered environment, stomatal conductance can be expected to be high, resulting in increased . A higher instantaneous elevates the overall Δ signature. Conversely, as soil moisture becomes scarce, stomatal conductance tends to decrease, lowering and, ultimately, the Δ signature . As such, the average Δ of all genotypes from a particular environment can provide an indication as to the severity of terminal drought. Unstressed plants typically register Δ values between 19 and 22 per mil . Being that the highest average Δ in this study was 18.2 per mil, there is evidence to suggest that each environment experienced a certain degree of terminal drought. On the other end of the spectrum, a severely stressed plant may register Δ values of 12 to 14 per mil . The lowest average Δ signature observed in the current study was 16.5 per mil, indicating high, but not severe, stress. In terms of values, the least severe terminal drought environment experienced average internal carbon concentrations of roughly 0.58, whereas the most severe saw average values of around 0.51 . A difference in average instantaneous of 7% over the period of grain carbon assimilation has the potential to result in large cumulative differences in biomass development. On the continuum from severe to moderate terminal drought stress, the environments included in this study are ordered as follows: Lind 2017, Pullman 2015, Pullman 2016, Pullman 2017, and Pendleton 2017.
Given that each environment experienced some level of terminal drought, it is likely that the Δ phenotype data was impacted by stomatal conductance, photosynthetic capacity, flowering time, root depth, and the translocation of stem reserves to varying degrees by genotype. The significant SNPs identified in this study may, therefore, also be associated with these contributing traits, excluding flowering time, as the addition of heading date as a covariate in the GWAS model helped to eliminate any associations due to genotypic differences in flowering time. Given the polygenic nature of each contributing trait [8, 47, 48], the dispersed genetic architecture for grain Δ revealed in this study is sensible. In a biparental mapping study involving grain Δ analysis on genotypes exposed to controlled terminal drought conditions, Peleg et al.  detected QTL on chromosomes 1A, 4A, 5A, 5B, 6B, and 7B. Both of the regions on 5B explained 27% of the observed variance. In the current study, chromosome 5B also appeared to be influential to Δ expression, with the greatest number of significant MTAs, as well as three repeated regions across environments. Notably, however, the phenotypic variance explained by these MTAs was much lower, the highest being just 7.5%.
While regions on chromosome 5B and others were found to repeat across environments, in general, the associated regions appeared to be largely environment specific. Such genotype-by-environment interactions also appeared to be evident in the grain Δ association analysis performed by Mora et al. , reporting all SNPs to be particular to each environment. Conversely, Rebetzke et al.  described several repeated QTL both across populations and between environments. However, none of the stable chromosomes (2B, 6D, 7A, and 7B) indicated were found to repeat across environments in this study.
Incidentally, phenotypic differences driven by genotype-specific reactions to the drought stress may also help to explain the relatively low heritability for Δ observed in this study. Broad-sense heritability for grain Δ is often relatively high, ranging from 0.40  to 0.83 . Intermediate values have been reported [26, 50]. Wu et al.  estimated a low broad-sense heritability for grain Δ of 0.23. The low heritability reported in this study is likely the result of high environmental variance, as evidenced by the significant differences in average Δ between environments, and the associated genotype-by-environment interactions. To improve heritability, it may be advisable to measure Δ on less stressed leaf tissue prior to the onset of terminal drought, as proposed by Condon et al. .
Grain Δ was determined to be a largely ineffective selection criterion for crop improvement in Washington target breeding environments. The low heritability combined with an inconsistent and generally weak correlation with grain yield makes selection for grain Δ unfavorable. The variation in environments greatly impacted the ability to use Δ as a selection tool. In areas with more consistent environmental factors, selection may be useful; however, in production regions with high variation between locations and years, Δ does not prove to be a valuable selection criterion. By measuring Δ on less stressed leaf tissue, it is possible that the genotype-by-environment interactions associated with grain Δ may diminish, ultimately improving heritability. The genetic architecture of grain Δ was found to be characterized by multiple, small effect, and largely environment-specific MTAs, further complicating the use of grain Δ as a selection criterion. These results suggest that grain Δ from wheat exposed to terminal drought is not amenable to MAS, as no single marker resulted in a large phenotypic response.
Jayfred V. Godoy’s present address is InterGrain Inc., Bibra Lake, WA, Australia 6163.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.
LSD analyzed the data and drafted the manuscript; JVG edited the manuscript; AHC edited the manuscript and obtained funding for the project.
This project was supported by the National Institute of Food and Agriculture, U.S. Department of Agriculture, under award 2016-68004-24770 and Hatch Project 1014919.
Supplementary 1. Table S1: monthly and total precipitation (mm) over the wheat growing season for five environments including Pullman (2015, 2016, and 2017), Lind 2017, and Pendleton 2017.
Supplementary 2. Table S2: variance components for calculation of broad-sense heritability () including total phenotypic variance (), environmental variance (), and genetic variance () for grain yield (GY, t ha-1), Δ, plant height (HT, cm), and heading date (HD, Julian) of a winter wheat panel grown in five Pacific Northwest environments.
Supplementary 3. Table S3: markers associated with carbon isotope discrimination (Δ) repeated in more than one environment.
Supplementary 4. Table S4: determination of linkage blocks by pairwise evaluation of markers associated with carbon isotope discrimination (Δ) that are equal to or less than 4 cM apart and are in significant linkage disequilibrium ().
Supplementary 5. Table S5: markers associated with grain yield from a panel of 480 advanced soft white winter wheat lines adapted to the Pacific Northwest and grown in five environments.
Supplementary 6. Figure S1: quantile-quantile plots for association analysis of carbon isotope discrimination (Δ) including principle components 1 and 2, and covariates for plant height and heading date provided a significant () Pearson correlation with Δ evaluated for each location-year including (A) Pullman 2015, (B) Pullman 2016, (C) Pullman 2017, (D) Lind 2017, (E) Pendleton 2017, and (F) best linear unbiased prediction (BLUP).
Supplementary 7. Figure S2: quantile-quantile plots for association analysis of grain yield including no covariates for each location-year including (A) Pullman 2015, (B) Pullman 2016, (C) Pullman 2017, (D) Lind 2017, (E) Pendleton 2017, and (F) best linear unbiased prediction (BLUP).
- H. Budak, M. Kantar, and K. Y. Kurtoglu, “Drought tolerance in modern and wild wheat,” The Scientific World Journal, vol. 2013, Article ID 548246, 16 pages, 2013.
- W. H. Pfeiffer, R. M. Trethowan, M. van Ginkel, M. I. Ortiz, and S. Rajaram, “Breeding for abiotic stress tolerance in wheat,” in Abiotic Stresses: Plant Resistance through Breeding and Molecular Approaches, M. Ashraf and P. J. C. Harris, Eds., pp. 401–489, Haworth Press, New York, NY, USA, 2005.
- R. T. Fortenberry and T. P. Nadreau, “Contribution of wheat production to the Washington economy,” Tech. Rep., IMPACT Center, Washington State University, 2017.
- M. M. Fontaine and A. C. Steinemann, “Assessing vulnerability to natural hazards: impact-based method and application to drought in Washington State,” Natural Hazards Review, vol. 10, no. 1, pp. 11–18, 2009.
- K. McLain and J. Hancock, Interim Report: 2015 Drought and Agriculture, Washington State Department of Agriculture, 2015.
- P. A. Knapp, P. T. Soulé, and H. D. Grissino-Mayer, “Occurrence of sustained droughts in the interior Pacific Northwest (A.D. 1733–1980) inferred from tree-ring data,” Journal of Climate, vol. 17, no. 1, pp. 140–150, 2004.
- S. Rauf, J. M. Al-Khayri, M. Zaherieva, P. Monneveux, and F. Khalil, “Breeding strategies to enhance drought tolerance in crops,” in Advances in Plant Breeding Strategies: Agronomic, Abiotic and Biotic Stress Traits, J. Al-Khayri, S. Jain, and D. Johnson, Eds., pp. 397–445, Springer, Cham, 2015.
- P. Gupta, H. S. Balyan, and V. Gahlaut, “QTL analysis for drought tolerance in wheat: present status and future possibilities,” Agronomy, vol. 7, no. 1, p. 5, 2017.
- R. A. Richards, G. J. Rebetzke, A. G. Condon, and A. F. van Herwaarden, “Breeding opportunities for increasing the efficiency of water use and crop yield in temperate cereals,” Crop Science, vol. 42, no. 1, p. 111, 2002.
- A. G. Condon, R. A. Richards, G. J. Rebetzke, and G. D. Farquhar, “Improving intrinsic water-use efficiency and crop yield,” Crop Science, vol. 42, no. 1, 2002.
- A. G. Condon, G. D. Farquhar, G. J. Rebetzke, and R. A. Richards, “The application of carbon isotope discrimination in cereal improvement for water-limited environments,” in Drought Adaptation in Cereals, pp. 171–219, Food Product Press, 2006.
- G. D. Farquhar and R. A. Richards, “Isotopic composition of plant carbon correlates with water-use efficiency of wheat genotypes,” Australian Journal of Plant Physiology, vol. 11, no. 6, pp. 539–552, 1984.
- G. D. Farquhar, J. R. Ehleringer, and K. T. Hubick, “Carbon isotope discrimination and photosynthesis,” Annual Review of Plant Physiology and Plant Molecular Biology, vol. 40, no. 1, pp. 503–537, 1989.
- E. Brugnoli and G. D. Farquhar, “Photosynthetic fractionation of carbon isotopes,” in Photosynthesis, R. C. Leegood, T. D. Sharkey, and S. Caemmerer, Eds., vol. 9 of Advances in Photosynthesis and Respiration, pp. 399–434, Springer, Dordrecht, 2000.
- A. G. Condon and R. A. Richards, “28 - Exploiting genetic variation in transpiration efficiency in wheat: an agronomic view,” in Stable Isotopes and Plant Carbon-Water Relations, pp. 435–450, Academic Press, 1993.
- G. J. Rebetzke, A. G. Condon, R. A. Richards, and G. D. Farquhar, “Selection for reduced carbon isotope discrimination increases aerial biomass and grain yield of rainfed bread wheat,” Crop Science, vol. 42, no. 3, 2002.
- A. Korte and A. Farlow, “The advantages and limitations of trait analysis with GWAS: a review,” Plant Methods, vol. 9, no. 1, article 29, 2013.
- L. T. Burghardt, N. D. Young, and P. Tiffin, “A guide to genome-wide association mapping in plants,” Current Protocols in Plant Biology, vol. 2, pp. 22–38, 2017.
- B. J. Soto-Cerda and S. Cloutier, “Association mapping in plant genomes,” in Genetic Diversity in Plants, IntechOpen, 2012.
- B. D. Singh and A. K. Singh, Marker-Assisted Plant Breeding: Principles and Practices, Springer, 2015.
- R. Bernardo, “Molecular markers and selection for complex traits in plants: learning from the last 20 years,” Crop Science, vol. 48, no. 5, pp. 1649–1664, 2008.
- Z. Peleg, T. Fahima, T. Krugman et al., “Genomic dissection of drought resistance in durum wheat × wild emmer wheat recombinant inbreed line population,” Plant, Cell & Environment, vol. 32, no. 7, pp. 758–779, 2009.
- G. J. Rebetzke, A. G. Condon, G. D. Farquhar, R. Appels, and R. A. Richards, “Quantitative trait loci for carbon isotope discrimination are repeatable across environments and wheat mapping populations,” Theoretical and Applied Genetics, vol. 118, no. 1, pp. 123–137, 2008.
- J. L. Araus, D. Villegas, N. Aparicio et al., “Environmental factors determining carbon isotope discrimination and yield in durum wheat under Mediterranean conditions,” Crop Science, vol. 43, no. 1, pp. 170–180, 2003.
- P. K. Gupta, H. S. Balyan, V. Gahlaut, and P. L. Kulwal, “Phenotyping, genetic dissection, and breeding for drought and heat tolerance in common wheat: status and prospects,” in Plant Breeding Reviews, J. Janick, Ed., pp. 85–168, John Wiley and Sons, 2012.
- G. Zhang, R. Aiken, and T. J. Martin, “Relationship between carbon isotope discrimination and grain yield of rainfed winter wheat in a semi-arid region,” Euphytica, vol. 204, no. 1, pp. 39–48, 2014.
- A. Pozo, A. Yáñez, I. A. Matus et al., “Physiological traits associated with wheat yield potential and performance under water-stress in a Mediterranean environment,” Frontiers in Plant Science, vol. 7, p. 987, 2016.
- P. Monneveux, M. P. Reynolds, R. Trethowan, H. Gonzalez-Santoyo, R. J. Pena, and F. Zapata, “Relationship between grain yield and carbon isotope discrimination in bread wheat under four water regimes,” European Journal of Agronomy, vol. 22, no. 2, pp. 231–242, 2005.
- A. G. Condon and A. E. Hall, “3 - Adaptation to diverse environments: variation in water-use efficiency within crop species,” in Ecology in Agriculture, pp. 79–116, Academic Press, 1997.
- W. F. Schillinger, R. I. Papendick, S. O. Guy, P. E. Rasmussen, and C. van Kessel, “Dryland cropping in the western United States,” in Dryland Agriculture, G. A. Peterson, P. W. Unger, and W. A. Payne, Eds., pp. 365–393, Agronomy Monograph no 23. ASA, CSSA, and SSSA, Madison, WI, USA, 2nd edition, 2006.
- P. S. Froese and A. H. Carter, “Single nucleotide polymorphisms in the wheat genome associated with tolerance of acidic soils and aluminum toxicity,” Crop Science, vol. 56, no. 4, pp. 1662–1677, 2016.
- R. E. Allan, C. J. J. Peterson, G. L. Rubenthaler, R. F. Line, and D. E. Roberts, “Registration of ‘Madsen’ wheat,” Crop Science, vol. 29, no. 6, 1989.
- A. E. Hall, R. A. Richards, A. G. Condon, G. C. Wright, and G. D. Farquhar, “Carbon isotope discrimination and plant breeding,” in Plant Breeding Reviews, J. Janick, Ed., pp. 81–113, John Wiley and Sons, 1994.
- G. Acquaah, Principles of Plant Genetics and Breeding, Wiley-Blackwell, 2nd edition, 2012.
- F. Breseghello and M. Sorrells, “Association mapping of kernel size and milling quality in wheat (Triticum aestivum L.) cultivars,” Genetics, vol. 172, no. 2, pp. 1165–1177, 2006.
- K. L. Jernigan, J. V. Godoy, M. Huang et al., “Genetic dissection of end-use quality traits in adapted soft white winter wheat,” Frontiers in Plant Science, vol. 9, 271 pages, 2018.
- P. J. Bradbury, Z. Zhang, D. E. Kroon, T. M. Casstevens, Y. Ramdoss, and E. S. Buckler, “TASSEL: software for association mapping of complex traits in diverse samples,” Bioinformatics, vol. 23, no. 19, pp. 2633–2635, 2007.
- Y. Tang, X. Liu, J. Wang et al., “GAPIT version 2: an enhanced integrated tool for genomic association and prediction,” The Plant Genome, vol. 9, no. 2, 2016.
- J. M. Bland and D. G. Altman, “Multiple significance tests: the Bonferroni method,” BMJ, vol. 310, no. 6973, p. 170, 1995.
- W. S. Bush and J. H. Moore, “Chapter 11: genome-wide association studies,” PLoS Computational Biology, vol. 8, no. 12, article e1002822, 2012.
- S. Wang, D. Wong, K. Forrest et al., “Characterization of polyploid wheat genomic diversity using a high-density 90,000 SNP array,” Plant Biotechnology Journal, vol. 12, no. 6, pp. 787–796, 2014.
- J. L. Araus, T. Amaro, J. Casadesu’s, A. Asbati, and M. M. Nachit, “Relationships between ash content, carbon isotope discrimination and yield in durum wheat,” Australian Journal of Plant Physiology, vol. 25, no. 7, pp. 835–842, 1998.
- S. C. Misra, R. Randive, V. S. Rao, M. S. Sheshshayee, R. Serraj, and P. Monneveux, “Relationship between carbon isotope discrimination, ash content and grain yield in wheat in the peninsular zone of India,” Journal of Agronomy and Crop Science, vol. 192, no. 5, pp. 352–362, 2006.
- O. Merah, E. Deléens, I. Souyris, M. M. Nachit, and P. Monneveux, “Stability of carbon isotope discrimination and grain yield in durum wheat,” Crop Science, vol. 41, no. 3, pp. 677–681, 2001.
- X. Xu, H. Yuan, S. Li, R. Trethowan, and P. Monneveux, “Relationship between carbon isotope discrimination and grain yield in spring wheat cultivated under different water regimes,” Journal of Integrative Plant Biology, vol. 49, no. 10, pp. 1497–1507, 2007.
- F. Mora, D. Castillo, B. Lado et al., “Genome-wide association mapping of agronomic traits and carbon isotope discrimination in a worldwide germplasm collection of spring wheat using SNP markers,” Molecular Breeding, vol. 35, no. 2, p. 69, 2015.
- Y. Dong, J. Liu, Y. Zhang et al., “Genome-wide association of stem water soluble carbohydrates in bread wheat,” PLoS One, vol. 11, no. 11, article e0164293, 2016.
- F. Shahinnia, J. Le Roy, B. Laborde et al., “Genetic association of stomatal traits and yield in wheat grown in low rainfall environments,” BMC Plant Biology, vol. 16, no. 1, article 150, 2016.
- A. G. Condon and R. A. Richards, “Broad sense heritability and genotype x environment interaction for carbon isotope discrimination in field-grown wheat,” Australian Journal of Agricultural Research, vol. 43, no. 5, pp. 921–934, 1992.
- B. Ehdaie and J. G. Waines, “Genetic analysis of carbon isotope discrimination and agronomic characters in a bread wheat cross,” Theoretical and Applied Genetics, vol. 88, no. 8, pp. 1023–1028, 1994.
- X. Wu, X. Chang, and R. Jing, “Genetic analysis of carbon isotope discrimination and its relation to yield in a wheat doubled haploid population,” Journal of Integrative Plant Biology, vol. 53, pp. 719–730, 2011.