Applied phenomics and genomics for improving barley yellow dwarf resistance in winter wheat

Abstract

Barley yellow dwarf is one of the major viral diseases of cereals. Phenotyping barley yellow dwarf in wheat is extremely challenging due to similarities to other biotic and abiotic stresses. Breeding for resistance is additionally challenging as the wheat primary germplasm pool lacks genetic resistance, with most of the few resistance genes named to date originating from a wild relative species. The objectives of this study were to (1) evaluate the use of high-throughput phenotyping to improve barley yellow dwarf assessment; (2) identify genomic regions associated with barley yellow dwarf resistance; and (3) evaluate the ability of genomic selection models to predict barley yellow dwarf resistance. Up to 107 wheat lines were phenotyped during each of 5 field seasons under both insecticide treated and untreated plots. Across all seasons, barley yellow dwarf severity was lower within the insecticide treatment along with increased plant height and grain yield compared with untreated entries. Only 9.2% of the lines were positive for the presence of the translocated segment carrying the resistance gene Bdv2. Despite the low frequency, this region was identified through association mapping. Furthermore, we mapped a potentially novel genomic region for barley yellow dwarf resistance on chromosome 5AS. Given the variable heritability of the trait (0.211–0.806), we obtained a predictive ability for barley yellow dwarf severity ranging between 0.06 and 0.26. Including the presence or absence of Bdv2 as a covariate in the genomic selection models had a large effect for predicting barley yellow dwarf but almost no effect for other observed traits. This study was the first attempt to characterize barley yellow dwarf using field-high-throughput phenotyping and apply genomic selection to predict disease severity. These methods have the potential to improve barley yellow dwarf characterization, additionally identifying new sources of resistance will be crucial for delivering barley yellow dwarf resistant germplasm.

barley yellow dwarf (BYD), high-throughput phenotyping (HTP), Triticum aestivum, virus, resistance, tolerance, genome-wide association mapping (GWAS), genomic selection (GS)

Introduction

Wheat (Triticum aestivum L.) is one of the most essential food crops in the world and is constantly threatened by biotic stresses (Savary et al. 2019). Among the most important viral stresses is barley yellow dwarf (BYD). This disease is widespread across the world and transmitted by aphids (Shah et al. 2012), and can cause significant yield reductions in susceptible cultivars. In Kansas, BYD is the fourth most significant wheat disease in terms of average estimated yield losses with an average yield loss of approximately 1% estimated over the past 20 years (Hollandbeck et al. 2019), equivalent to a loss of more than $10 million per year. However, yield losses are highly variable ranging from 5% to 80% in a single field depending on the environment, management practices, the host, and the genetic background (Miller and Rasochová 1997; Perry et al. 2000; Gaunce and Bockus 2015). Moreover, the wide host range and the complex lifestyle of its vectors make BYD extremely difficult to manage, and different management strategies (e.g. planting date and control of vector populations) are inconsistent depending on climate and location (Bockus et al. 2016). Thus, in many production environments, particularly in the Central and Eastern regions of Kansas, BYD is often the most economically impactful disease.

BYD disease symptoms are highly variable depending on the crop, variety, time, and developmental stage when the infection occurs, aphid pressure, and environmental conditions (Shah et al. 2012; Choudhury et al. 2019b). BYD characterization in the field is extremely challenging as the symptoms can easily be confused with other viral disease symptoms such as wheat streak mosaic virus symptoms, nutrient deficiencies, or environmental stresses such as waterlogging (Shah et al. 2012). Typical BYD symptoms can be observed at all levels of plant organization—leaf, roots, and flowers. Leaf discoloration in shades of yellow, red, or purple, specifically starting at the tip of the leaf and spreading from the margins toward the base is common as well as a reduction in chlorophyll content (Jensen and Van Sambeek 1972; D’arcy 1995). Often the entire plant visually appears stunted or dwarfed from a reduction in biomass by reducing tiller numbers. Spike grain yield is decreased through a reduction in kernels per spike and kernel weight which also affects grain quality (Riedell et al. 2003; Choudhury et al. 2019b). Quality can be further reduced by a reduction in starch content (Peiris et al. 2019). Below ground effects of BYD have also been reported including reduced root growth (Riedell et al. 2003).

Currently, there is no simple solution to control BYD (Walls et al. 2019), however, the use of genetic resistance and tolerance is the most appealing and cost-effective option to control this disease (Comeau and Haber 2002; Choudhury et al. 2017, 2019b). Resistance and tolerance could be different genetic mechanisms, namely stopping virus replication and minimizing disease symptoms, respectively, but within this paper all mention of resistance includes both genetic resistance and tolerance. Breeding strategies involving genetic resistance can target either the aphids or the virus. Resistance to aphids can be achieved by 3 different strategies, antixenosis, antibiosis, or tolerance (Girvin et al. 2017). To date, most breeding efforts have been directed to the identification of viral tolerance, also known as “field resistance,” that refers to the ability of the plant to maintain yield under BYD infection and is associated with a reduction of symptoms of infection independent of the virus titer (Foresman et al. 2016). Field resistance has been reported to be polygenic, falling under the quantitative resistance class, where several genes with very small effects control the resistance response (Qualset et al. 1973, Cisar et al. 1982; Ayala et al. 2002; Choudhury et al. 2019a, 2019c).

Presently, no major gene conferring immunity or a strong resistant phenotype to BYD has been identified in bread wheat, and only 4 resistance genes have been described for BYD. Located on chromosome 7DS, Bdv1 is the only gene described from the primary pool of wheat and was originally identified in the wheat cultivar “Anza” (Singh et al. 1993). This gene provides resistance to some but not all the viruses that cause BYD (Ayala-Navarrete and Larkin 2011). The other 3 named genes were all introduced into wheat through wide crossing from intermediate wheatgrass (Thinopyrum intermedium) (Ayala et al. 2001; Zhang et al. 2009). Bdv2 and Bdv3 are both located on a translocation segment on wheat chromosome 7DL (Brettell et al. 1988; Sharma et al. 1995), while Bdv4 is located on a translocation segment on chromosome 2D (Larkin et al. 1995; Lin et al. 2007). Bdv2 was the first gene successfully introgressed in wheat breeding programs from the tertiary gene pool for BYD resistance (Banks et al. 1995) and deployed into varieties.

In addition to the 4 known resistance genes, other genomic regions associated with BYD resistance have been identified through genetic mapping. These regions have been described on nearly all wheat chromosomes but have not been genetically characterized (Ayala et al. 2002; Jarošová et al. 2016; Choudhury et al. 2019a, 2019b, 2019c). Of the described regions, most explain a minor proportion of the genetic variation (<15%) (Ayala et al. 2002; Choudhury et al. 2019a, 2019c) in biparental populations suggesting a potential upwardly biased estimate due to the Beavis effect (Xu 2003). Moreover, 2 recent studies have reported that some of these new genomic regions display additive effects (Choudhury et al. 2019a, 2019b). Additive genetic effects had already been reported in lines combining Bdv2 and Bdv4 (Jahier et al. 2009).

Taken together, research indicates that resistance genes to BYD in wheat are rare. With a lack of major genes and difficulty to characterize resistance in the wheat pool likely due to the polygenic nature of many small effect loci, identifying resistance has been limited. Nevertheless, breeding programs have devoted large efforts for breeding BYD resistance due to the economic importance of this disease, with some of the greatest success coming from wide crosses to the tertiary gene pool.

Breeding for BYD resistance can be improved by applying strategies for more effective evaluation and utilization of the identified resistance. To get a better understanding of BYD and its quantitative nature, consistent, and high-throughput methods are needed for the identification of resistant wheat lines for large-scale selection in breeding programs (Aradottir and Crespo-Herrera 2021). Effective selection of quantitative resistance with low heritability can be aided by high-throughput genotyping, high-throughput phenotyping (HTP), or a combination of both.

Access to high-density genetic markers at a very low-cost, owing to the rapid developments in DNA sequencing, have enabled breeding programs to apply molecular breeding for quantitative traits. Genomic selection (GS) is a powerful tool to breed for quantitative traits with complex genetic architecture and low heritability (e.g. yield, quality, and diseases such as Fusarium head blight), because it has greater power to capture loci with small effect compared with other marker-assisted selection strategies (Meuwissen et al. 2001; Poland and Rutkoski 2016). In addition to molecular data, HTP using unmanned aerial systems (UAS) or ground-based sensors is providing high density phenotypic data that can be incorporated into breeding programs to increase genetic gain (Haghighattalab et al. 2016; Crain et al. 2018; Wang et al. 2020). Using precision phenotyping for disease scoring can improve the capacity for rapid and nonbiased evaluation of large field-scale numbers of entries (Poland and Nelson 2011). Taken together, improvements in genomics and phenomics have the potential to aid breeding progress for BYD resistance.

In an effort to accelerate the development of resistant lines, we combined high throughput genotyping and phenotyping to assess BYD severity in a large panel of elite wheat lines. We evaluated the potential of HTP data to accurately assess BYD severity as well as identify genetic regions associated with BYD resistance and inform whole genome prediction to identify resistant lines.

Materials and methods

Plant material

A total of 381 different wheat genotypes were characterized for BYD resistance, including 30 wheat cultivars and 351 advanced breeding lines in field nurseries over 5 years (Supplementary Table 1). In each nursery, an unbalanced set of 52–107 wheat entries were evaluated including both cultivars and breeding lines (Table 1). The BYD susceptible cultivar “Art” and BYD resistant cultivar “Everest” were included in all the nurseries (seasons) as checks, and no other wheat genotype was in common between different seasons of evaluation.

Table 1.

Field experimental details for the 5 wheat nurseries.

Season	2015–2016	2016–2017	2017–2018	2018–2019	2019–2020
Location	Rocky Ford farm		Ashland Bottoms farm
	39°13′45.60″N, 96°34′41.21″W		39°07′53.76″N, 96°37′05.20″W
Planting date	2015 September 17	2016 September 12	2017 September 19	2018 September 17	2019 September 17
Number of entries	68	52	81	81	107
Number of plots	504	360	400	392	476
Field design	Split-plot with insecticide treatment as main factor effect and wheat genotype as secondary factor
Replications	3	3	2	2	2
Plot size	6 rows plots—1.5 m × 2.4 m
BYD evaluation	2016 April 28	2017 May 12	2018 May 19	2019 May 13	2020 May 19
Harvesting date	2016 June 20	2017 June 19	2018 June 23	2019 June 28	2020 June 25

Field experiments

Nurseries for BYD field-screening were conducted during 5 consecutive wheat seasons (2015–2016 to 2019–2020) (Table 1). Seasons 2015–2016 and 2016–2017 were conducted at Kansas State University (KSU) Rocky Ford experimental station (39°13′45.60″N, 96°34′41.21″W), while the 2017–2018, 2018–2019, and 2019–2020 nurseries were planted at KSU Ashland Bottoms experimental station (39°07′53.76″N, 96°37′05.20″W). The nurseries were established for natural BYD infections by planting in mid-September, about 3 weeks earlier than the normal planting window. The susceptible cultivar “Art” was planted as a spreader plot in the borders and as a control check plot also with the resistant cultivar “Everest.” The experimental unit was 1.5 m × 2.4 m with a 6-row plot on 20 cm row spacing.

A split-plot field design with 2 or 3 replications was used where the main plot was insecticide treatment, and the split plot was the wheat genotype. Three replications were used for proof of concept during the first 2 seasons but then 2 replications were chosen as a balance of space and number of entries for the following seasons. For the treated replications the seed was treated at planting with Gaucho XT (combination of insecticide and fungicide) at a rate of 0.22 ml/100 g of seed, followed with foliar insecticide applications starting approximately 2–3 weeks after planting through heading. Depending on field conditions, spray treatments were conducted every 14–21 days if average air temperatures remained above 10°C. Foliar insecticides were applied to the treated replications in a spray volume of 280.5 l/ha using a Bowman MudMaster plot sprayer equipped with TeeJet Turbo TwinJet tips. Insecticide applications consisted of a rotation of Warrior II, Lorsban, and Mustang Max at rates of 0.14, 1.17, and 0.29 l/ha, respectively. For the control insecticide treatment (untreated), the seed was treated with Raxil MD (fungicide) at a rate of 0.28 ml/100 g of seed, and no foliar insecticide applications were applied. Foliar fungicide Nexicor was applied to the whole experiment at a rate of 0.73 l/ha, at both planting and heading, to control all other diseases so the main disease pressure was focused on BYD.

Phenotypic data

Individual plots were assessed for (1) BYD severity characterized as the typical visual symptoms of yellowing or purpling on leaves using a 0–100% visual scale, determined directly after spike emergence by recording the proportion of the plot exhibiting the symptoms (Table 1); (2) manual plant height (PTHT_M, m); and (3) grain yield (GY, tons/ha). Experimental plots were harvested using a Kincaid 8XP plot combine (Kincaid Manufacturing., Haven, KS, USA). Grain weight, grain moisture and test weight measurements for each plot was recorded using a Harvest Master Classic GrainGage and Mirus harvest software (Juniper Systems, Logan, UT, USA). Visual phenotypic assessment was recorded using the Field Book phenoapp (Rife and Poland 2014).

High-throughput phenotyping

To compliment the manually recorded phenotypic data, we applied HTP using a ground-based proximal sensing platform or an UAS (Table 2). Seasons 2015–2016 and 2016–2017 were characterized by the ground platform as described in Barker et al. (2016) and Wang et al. (2018). For the other 3 seasons, we used a quadcopter DJI Matrice 100 (DJI, Shenzhen, China) carrying a MicaSense RedEdge-M multispectral camera (MicaSense Inc., USA). The HTP data were collected on multiple dates throughout the growth cycle from stem elongation to ripening (GS 30–90; Zadoks et al. 1974) (Table 2). Flight plans were created using CSIRO mission planner application and missions were executed using the Litchi Mobile App (VC Technology Ltd., UK; https://uavmissionplanner.netlify.app/) for DJI Matrice100. The aerial image overlap rate between 2 geospatially adjacent images was set to 80% both sequentially and laterally to ensure optimal orthomosaic photo stitching quality. All UAS flights were set at 20 m above ground level at 2 m/s and conducted within 2 h of solar noon. To improve the geospatial accuracy of orthomosaic images, white square tiles with a dimension of 0.30 m × 0.30 m were used as ground control points and were uniformly distributed in the field experiment before image acquisition and surveyed to centimeter-level resolution using the Emlid REACH RS+ Real-Time Kinematic Global Navigation Satellite System unit (Emlid Ltd, Hong Kong, China).

Table 2.

Dates of high-throughput phenotypic data collection and details of image acquisition in the 5 wheat nurseries screened for BYD, Kansas, USA (2015–2020).

Season	2015–2016	2016–2017	2017–2018	2018–2019	2019–2020
UAS platform	PheMU		DJI Matrice 100
Imaging sensor	Multiple digital single-lens reflex (DSLR) cameras		MicaSense RedEdge-M
Flight/pass speed	0.3–0.5 m/s		2 m/s
Flight dates	2016-03-31 2016-04-07 2016-04-14 2016-05-06	2017-03-28 2017-04-13 2017-05-01 2017-05-09 2017-05-21 2017-05-23 2017-05-30 2017-06-05 2017-06-13	2018-03-30 2018-04-04 2018-04-12 2018-04-19 2018-04-23 2018-05-16 2018-06-13	2019-04-01 2019-04-09 2019-04-19 2019-04-26 2019-05-02 2019-05-10 2019-05-15 2019-05-23 2019-05-31 2019-06-05 2019-06-12 2019-06-17	2020-03-20 2020-04-11 2020-04-23 2020-05-03 2020-05-19 2020-06-05 2020-06-11
Flight/pass altitude	0.5 m above the canopy		20 m AGL
In-air flight duration	NA		∼11–14 min

An automated image processing pipeline (Wang et al. 2020) was used to generate the orthomosaics and extract plot-level plant height [PTHT_D (m); Singh et al. 2019] and normalized difference vegetation index (NDVI) (Rouse et al. 1974), calculated as:

$$\mathrm{NDVI}=\mathrm{ }\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}},$$

(1)

where NIR and Red are the near-infrared and red bands of the multispectral images and NDVI is the output image. Both traits were selected based on potential BYD characterization where the most typical BYD symptoms include chlorosis and stunting of the plants, thus, influencing NDVI and PTHT.

Statistical data analyses

First, the adjusted mean best linear unbiased estimator (BLUE) was calculated for each entry for all the different traits for each season (Supplementary Table 1), using the following model:

$$y_{\mathit{ijklm}}={\mu +G}_i+T_j+{GT}_{ij}+R_{k\left(j\right)}+B_{l\left(kj\right)}+C_{m\left(kj\right)}+e_{\mathit{ijklm}},$$

(2)

where $y_{\mathit{ijklm}}$ is the phenotype for the trait of interest, $\mu$ is the overall mean, $G_i\mathrm{ }$is the fixed effect of the $i\mathrm{th}$ entry (genotype), $T_j\mathrm{ }$is the fixed effect of the $j\mathrm{th}$ insecticide treatment, ${GT}_{ij}\mathrm{ }$is the fixed effect of the interaction between the $i\mathrm{th}$ entry and the $j\mathrm{th}$ insecticide treatment (genotype by treatment effect), $R_{k\left(j\right)}\mathrm{ }$is the random effect of the $k\mathrm{th}$ replication nested within the $j\mathrm{th}$ insecticide treatment and distributed as iid $R_{k\left(j\right)}\sim N(0,{\sigma }_R^2)$⁠, $B_{l\left(kj\right)}\mathrm{ }$is the random effect of the $l\mathrm{th}$ row nested within the $k\mathrm{th}$ replication and $j\mathrm{th}$ treatment distributed as iid $B_{l\left(kj\right)}\sim N(0,{\sigma }_B^2)$⁠, $C_{m\left(kj\right)}\mathrm{ }$is the random effect of the $m\mathrm{th}\mathrm{ }$column nested within the $k\mathrm{th}$ replication and $j\mathrm{th}$ treatment and assumed distributed as iid $C_{m\left(kj\right)}\sim N(0,{\sigma }_C^2)$⁠, and $e_{\mathit{ijklm}}$ is the residual for the $\mathit{ijklm}\mathrm{th}$ plot and distributed as iid $e_{\mathit{ijklm}}\sim N(0,{\sigma }_e^2)$⁠. The “lme4” R package (Bates et al. 2015) was used for fitting the models.

The BLUEs were used to inspect trait distributions and to calculate Pearson’s correlations between all traits. In addition, BLUE values were used to calculate the reduction in GY for each entry as the difference of GY between the untreated and insecticide treated main plots. This variable reflects the level of BYD resistance of each entry, and it was used to perform GWAS and GS analyses.

For NDVI and PTHT_D, the plot-level observed values extracted for the different phenotypic dates were fitted to a logistic nonlinear regression model (Fox and Weisberg 2011) as,

$$y=\mathrm{ }\frac{{\theta }_1}{1+\mathrm{ }{e}^{-({\theta }_2+{\theta }_{3}x)}}+\mathrm{ }\epsilon ,$$

(3)

where $y$ is the phenotype for the trait of interest at the time-point $x$ measured as days after January 1, ${\theta }_1$ is the maximum value (upper asymptote) represented by the final PTHT or maximum achieved NDVI, ${\theta }_2$ is the inflection point that represents the greatest rate of change in the growth curve, either senescence for NDVI or height of growth, ${\theta }_3$ is the lag phase or onset of senescence or growth rate from time $x$ where x is the calendar day of the year since January 1, and $\epsilon$ is the residual error (Supplementary Fig. 1). The “nlme” R package was used for model fitting (Pinheiro et al. 2015). The model parameters obtained for each trait (⁠${\theta }_{1N\mathit{DVI}}$⁠, ${\theta }_{2N\mathit{DVI}}$⁠, ${\theta }_{3N\mathit{DVI}}$⁠, ${\theta }_{1{\mathit{PTHT}}_D}$⁠, ${\theta }_{2{\mathit{PTHT}}_D}$⁠, and ${\theta }_{3{\mathit{PTHT}}_D}$⁠) were used in addition to the other phenotypic traits to calculate BLUEs, distributions, correlations, and BLUPs.

Secondly, we used a mixed linear model to calculate the best linear unbiased predictors (BLUPs) for each entry in each nursery (season) (Supplementary Table 1), using the same model as described in Equation (2) but defining $G_i$⁠, $T_j$⁠, and ${GT}_{ij}$ as random effects. BLUPs were used because of the unbalanced nature of the data. The BLUPs calculated for each season were then combined for GWAS and GS. Furthermore, we calculated broad-sense heritability on a line-mean basis by splitting the data based on whole plot treatment for insecticide treatments as:

$$H^2=\frac{{\sigma }_G^2}{{\sigma }_G^2+({\sigma }_e^2/\mathrm{r})},$$

(4)

where ${\sigma }_G^2$ is the genotypic variance, ${\sigma }_e^2$ is the residual error variance, and $r$ is the number of replications.

Genotypic data

A total of 346 wheat entries were genotyped using genotyping-by-sequencing (GBS) (Poland et al. 2012) and sequenced on an Illumina Hi Seq2000. Single nucleotide polymorphisms (SNPs) were called using Tassel GBSv2 pipeline (Glaubitz et al. 2014) and anchored to the Chinese Spring genome assembly v1.0 (IWGSC et al. 2018). SNP markers with minor allele frequency <0.01, missing data >85%, or heterozygosity >15% were removed from the analysis. After filtering, we retained 29,480 SNPs markers that were used to investigate the population structure through principal component analysis (PCA), genome-wide association analysis (GWAS), and GS. In addition, GBS data were used to run a bioinformatics pipeline to predict the presence or absence of the translocated segment on chromosome 7DL carrying the Bdv2 gene for each entry (Supplementary Table 1). The prediction was done based on a modified alien prediction pipeline (Gao et al. 2021). Briefly, alien or wheat specific tags were counted in the 7DL region and tabulated using a training set of cultivars or lines that are known to be Bdv2 positive and negative. A simple classification was done based on alien to wheat tag counts ratios.

Genome-wide association analysis

The GWAS analysis was performed with a mixed linear model according to Yu et al. (2006) implemented in the “GAPIT” R package (Lipka et al. 2012) that included the first 2 principal components to account for population structure as fixed effects and the individuals to explains familial relatedness as random effects,

$$y=\mathbf{X}\beta +\mathbf{Z}{u}_i+e,$$

(5)

where y is the vector of phenotypic BLUPs, $\mathbf{X}$ and $\mathbf{Z}$ are the incidence matrix of $\beta$ and $u_i$⁠, respectively, with $u_i$ assumed ∼ N (0, 2${\mathbf{K}}_i{\sigma }_i^2)$ where $\mathbf{K}$ is the individual kinship matrix, and $e$ is the vector of random residual effects with ∼ N (0, I${\sigma }_e^2)$⁠, where $\mathrm{I}$ is the identity matrix and ${\sigma }_e^2$ is the unknown residual variance. The false discovery rate correction with an experimental significance level value of 0.01 was used to assess marker-trait associations. Manhattan plots were generated with “CMplot” package in R software (Yin 2020). PCA using GBS-SNPs was performed in R language. Eigenvalues and eigenvectors were computed with “e” function using “A.mat” function and the “mean” imputation method of “rrBLUP” package (Endelman 2011). To declare a quantitative trait locus (QTL) we considered only the regions having several SNP markers in linkage disequilibrium, clearly showing a peak. We did not consider regions with a single SNP above the significant threshold as a QTL.

Genomic selection

Using data from the 5 seasons, GS models using the genomic BLUP (G-BLUP) were developed to assess predictive ability. A 5-fold cross-validation method was used to assess model accuracy where the data set was split into 5 sets based on season, with 4 seasons forming the training set and the fifth season serving as prediction set. This process was repeated until all seasons were predicted. Along with predicting all other seasons from each season, a model was evaluated with a leave-two-out cross-validation strategy. This strategy was used to get a better mix of years with and without disease incidence, where the training population consisted of 3 seasons, and the remaining 2 seasons were predicted from the combined training population. The GS model was fitted with the training population using “rrBLUP” kin.blup function (Endelman 2011), the GS model equation was,

$$y=\mathbf{W}g+\epsilon ,$$

(6)

where $y$ is a vector of phenotypic BLUPs,$\mathrm{ }\mathbf{W}$ is the design matrix of $g$⁠, $g$ is the vector of genotypic values ∼ N (0, $\mathbf{K}{\sigma }_g^2),$ and $\epsilon$ is the vector of residual errors (Endelman 2011). Predictive ability was assessed using Pearson’s correlation (⁠$r$⁠) between the predicted value (G-BLUP) and the BLUP for the respective phenotype. In addition, for both GS strategies we also tested the effect of adding the genotype of the Bdv2 loci as a fixed effect cofactor, using the model,

$$y=\mu +\mathbf{X}\beta +\mathbf{W}g+\epsilon ,$$

(7)

which combines parameters described in Equation (5) and $\mathrm{X}$ is the matrix (n × 1) of individual observation for presence or absence of Bdv2 and $\beta$ is the fixed effect for the Bdv2 measurements.

Results

Phenotypic data

We analyzed 5 years of BYD field-screening nurseries (seasons 2015–2016 to 2019–2020) characterizing a total of 381 wheat lines. The disease pressure and the expression of BYD associated symptoms varied each season, however, we were able to observe a significant effect of the insecticide treatment in all seasons (Fig. 1). Across all seasons, BYD symptoms were lower on the insecticide treated plots and both PTHT_M and GY increased compared with the nontreated control. Season 2016–2017 had the most conducive conditions for BYD screening, resulting in high average severity and a larger difference between mean values for the treated vs untreated plots for all the collected traits (Fig. 1). There was general consistency in order across all seasons with the susceptible check “Art” ranked among the highest in BYD severity (Supplementary Fig. 2).

Fig. 1.

Adjusted phenotypic values for the traits collected manually for 5 different field seasons (2015–2016 to 2019–2020). a–e): Barley yellow dwarf severity (%) characterized as the typical visual symptoms of yellowing/purpling on leaves using a 0–100% visual scale; f–i) manual plant height/stunting (PTHT_M) (m), note that the trait was not recorded for the 2015–2016 season; and j–n) grain yield (tons/ha). The dashed line represents the mean value for the trait in each treatment.

Phenotypic correlations between the traits showed a negative correlation between BYD and GY for all the seasons and a negative or no correlation between BYD and PTHT_M (Supplementary Fig. 3). The same correlation trends were observed under insecticide treated and untreated plots. Broad-sense heritability was moderate to high for all the traits, ranging between 0.21 and 0.79 for the insecticide treated plots and between 0.41 and 0.84 for the untreated plots. Across all traits, the untreated insecticide replications showed higher $H^2$values, with season 2016–2017 showing the highest values (Fig. 2).

Fig. 2.

Broad-sense heritability of wheat phenotypic traits collected manually, including visual barley yellow dwarf (BYD) score, plant height (PTHT_M), and grain yield (GY) during 5 different field seasons under 2 insecticide treatments.

For the HTP data collected (Table 2), we obtained 3 different parameters (⁠${\theta }_1$⁠, ${\theta }_2$⁠, and ${\theta }_3$⁠) for both PTHT_D and NDVI after fitting a logistic regression model using the data collected during the experiments (2015–2016 season data were not included due to lack of data quality) (Supplementary Fig. 1). Correlations between these parameters and the phenotypic traits collected manually were different for all the traits (Supplementary Fig. 3). For the insecticide untreated plots, BYD resulted in a negative correlation with ${\theta }_{2N\mathit{DVI}}$ and a positive correlation with ${\theta }_{3N\mathit{DVI}}$⁠, in most of the field seasons. We did not find a clear correlation pattern between BYD and${\mathrm{ }\mathit{PTHT}}_D$⁠. For PTHT_M we detected a positive correlation with ${\theta }_{1{\mathit{PTHT}}_D}$across all seasons, and for GY we observed a positive correlation with ${\theta }_{1N\mathit{DVI}}$ and ${\theta }_{2N\mathit{DVI}}$⁠, and a negative correlation with ${\theta }_{3N\mathit{DVI}}$ (Supplementary Fig. 3).

Prediction of Bdv2 resistance gene

We used GBS data to genotype the Bdv2 resistance gene located on a translocation segment from intermediate wheatgrass on chromosome 7DL of bread wheat. In total, 33 of the 346 wheat lines carried the Th. intermedium chromosomal translocation with Bdv2 (Supplementary Table 1). Interestingly, 28 of these Bdv2 lines belonged to the same breeding cycle, entering the advanced yield nursery stage of the KSU breeding program in the 2017–2018 season. Furthermore, only 7 pedigrees are represented within the 28 Bdv2 entries, meaning that these lines are highly related. The remaining 5 Bdv2 lines were distributed in 2015–2016 (n = 3), 2018–2019 (n = 1), and 2019–2020 (n = 1), and none of the lines from the season 2016–2017 had the presence of Bdv2 (Supplementary Table 1).

Population structure

We studied the population structure of 346 wheat lines using 29,480 GBS-derived SNP markers. The PCA did not reveal a strong pattern of population structure (Fig. 3). Moreover, the variation explained by the first 2 principal components (4% and 3%, respectively) also supports the hypothesis of minimal population structure within a single breeding program. We observed that most of the wheat cultivars released by KSU breeding program were located outside the cluster grouping all the breeding lines (Fig. 3a). Lines with the presence of Bdv2 clustered together (Fig. 3b), likely due to a related pedigree to the original source, and we did not identify any evident pattern for BYD severity associated with the population structure (Fig. 3c).

Fig. 3.

Scatterplot of the first two principal component axis, made from principal component analysis on the marker matrix, n = 357 wheat lines, markers = 29,480. Each data point represents an individual wheat line that is color-coded by (a) breeding status; (b) prediction of Bdv2 presence/absence; and (c) adjusted mean for BYD severity (BYD BLUE) scored visually. Total variance explained by each principal component (PC) is listed on the axis.

Genome-wide association analysis

To investigate the genetic architecture of BYD we performed GWAS analyses for all collected traits using the BLUP values for 346 lines and 29,480 SNP markers. The first two principal components from PCA and the kinship matrix were included in the mixed model to account for population structure and genetic relatedness. We found significant marker-trait associations for BYD severity on chromosomes 5AS, 7AL, and 7DL (Fig. 4a). The highest peak was observed on the proximal end of chromosome 7DL, located at 571–637 Mbp. To test the hypothesis that this association was explained by the resistance gene Bdv2 (located on chromosome 7DL), we investigated the haplotypes defined by the 16 SNP markers associated with BYD severity and were able to identify 2 haplotypes that exactly matched the presence or absence of Bdv2 (Fig. 4a). This same region was mapped using BYD severity and the presence or absence of Bdv2 as a fixed covariate (Fig. 4b). This analysis (Fig. 4b) also detected a peak on chromosome 7AL. Lastly, we explored the effect of Bdv2 on both BYD BLUEs and BLUPs, and we observed that the presence of Bdv2 had a positive effect in reducing the disease severity by approximately 10% (Fig. 5a). The significant peak on chromosome 5AS, located at 46–103 Mbp, was explained by 10 SNP markers, comprising 2 main haplotypes, one of them associated with reduced BYD severity (Fig. 5b). When we combined the different 5AS haplotypes with Bdv2, we observed that the presence of Bdv2 had a positive effect, reducing the levels of BYD when combined with both 5AS haplotypes (Fig. 5c), and suggesting an additive effect. In addition, there was no significant difference in BYD reduction between the 2 5AS haplotypes combined with presence of Bdv2 (Fig. 5d). Compared to the associations found for Bdv2 (Fig. 4b), we did not find any strong evidence of marker trait associations for the other evaluated traits (Supplementary Fig. 4).

Fig. 4.

Manhattan plots showing the marker-trait associations using 346 wheat accessions and 29,480 SNP markers obtained with genotyping-by-sequencing (GBS) for (a) BYD severity and (b) presence/absence of Bdv2 resistance gene. The 21 labeled wheat chromosomes with physical positions are on the x-axis and y-axis is the –log10 of the P-value for each SNP marker. Horizontal dashed lines represent the false discovery rate threshold at 0.01 level and highlighted data points above the threshold represent SNPs significantly associated with the trait. In (a), the length of the region and the haplotypes defined by the significant SNP markers is displayed.

Fig. 5.

Measurement of barley yellow dwarf disease severity in wheat based on certain haplotype effects were (a) the presence or absence of the translocation segment carrying the resistance gene Bdv2; (b) the 2 haplotypes for the significant region on chromosome 5AS; (c) the combination of 5A haplotypes with the presence or absence of Bdv2 gene; and (d) the 5A haplotypes combined with presence of Bdv2 resistant allele. Boxplots show the significant reduction of BYD disease severity by averaging the phenotypic best linear unbiased estimated (BLUE) values for the lines. Count is the number of wheat genotypes averaged in each group and mean is the mean BLUE value for the group.

Genomic selection

To evaluate the potential of GS to predict BYD disease severity, we fit several GS models to the phenotypic BLUPs of BYD, PTH_M, and reduction in GY. Across all traits, to determine predictive ability we used a 5-fold cross validation where prediction ability ranged from −0.08 to 0.26. There was relatively good predictive ability for BYD severity ranging between 0.06 and 0.26, in comparison with PTHT_M and reduction in GY resulting in a lower range from 0.02 to 017 and −0.08 to 0.2, respectively (Fig. 6). Evaluating the composition of the training population, we observed that when including 2016–2017 season, prediction abilities were the highest for BYD but the lowest for the other 2 traits, implying that season 2016–2017 was either a good season to train the prediction models or a difficult season to predict based on available data.

Fig. 6.

Genomic selection model predictive ability where each column represents one trait, and each row shows the conformation of the training population including size of training and testing population and number of lines with presence of Bdv2 resistance gene. The value in each cell represents the predictive ability which is the correlation between the GS predicted value (GBLUP) and the phenotypic BLUP.

To further investigate the ability of GS to predict BYD, we tested GS models using a leave-two-out strategy, where 2 seasons were excluded from the training population and used as the testing population. We fitted GS models for all possible 2-season combinations. This strategy resulted in slightly smaller training populations which decreased overall predictive ability (Fig. 6). This result was evident for BYD predictions where excluding 2 seasons had a larger negative impact.

Lastly, we evaluated the effect of adding information about the genotype of the Bdv2 resistance gene as a phenotypic fixed covariate into the GS models. There were differences in the effect of Bdv2 on the predictive ability across BYD severity, PTHT_M, and GY, showing a large effect for predicting BYD but almost no effect for PTHT_M and reduction in GY (Fig. 6). Including the presence or absence of Bdv2 as a covariate had a major effect on the predictive ability. For example, in the 2017–2018 season which had the highest proportion of lines with presence of Bdv2, excluding the covariate, resulted in a drop of the predictive ability from 0.09 to −0.21.

Discussion

Phenotypic data

The success of breeding for BYD resistance is highly impacted by the ability to precisely characterize breeding material and disease symptoms. Even though BYD is spread worldwide, its incidence in a given year depends on several factors such as aphid pressure, planting date, and environmental conditions (e.g. temperature, rainfall, frost, etc.). In this study, we evaluated winter wheat advanced breeding lines during 5 seasons implementing a rigorous field-testing approach, that ultimately enabled us to consistently have plots contrasting with BYD infection and uninfected or low incident plots. Moreover, by using large yield-size plots we were able to calculate the reduction in GY and use this parameter as an estimate of field resistance.

The expression of BYD symptoms, however, was highly inconsistent during the different seasons. Seasons 2015–2016 and 2016–2017 showed the best expression of the disease symptoms, supported by the wide range of BYD severity between treated and untreated replications (Fig. 1). Interestingly, both these seasons were conducted in the same experimental field (Table 1), suggesting that this location could favor the development of BYD. Moreover, weather conditions were variable for all the seasons, suggesting that these had a huge impact on the disease occurrence. While temperature records were similar for all the seasons, precipitation records did show some differences. Season 2017–2018 was dryer than normal, with 34% less precipitation than the 30 years historical average (1981–2010). On the other hand, season 2018–2019 was wetter than normal, with 58% more precipitation than the 30 years historical average (Supplementary Table 2).

High-throughput phenotyping

Evaluating BYD resistance using visual phenotypic selection can be challenging due to the complex nature of the disease and rater variability (Poland and Nelson 2011). The use of HTP with UAS is gaining popularity within plant breeding programs because it further improves selection intensity and accuracy compared with conventional phenotyping. Accurate phenotyping is crucial for understanding the genetic basis of quantitative and complex traits such as BYD severity. In this study, we used HTP to complement the visual BYD scoring. This phenotyping methodology improved our capacity for rapid, nondestructive, and nonbiased evaluation of large field-scale numbers of entries for BYD resistance. We were able to observe strong correlation patterns between visual BYD severity and HTP derived parameters (Supplementary Fig. 3); although, none of the HTP traits collected in this study had a common genetic base with BYD severity (Fig. 4; Supplementary Fig. 4). This uncertainty could be raised by 3 main technical reasons. Firstly, in this study a 5-band multispectral camera was used to capture spectral information reflected from canopies. Limitation due to the spectral resolution might restrict the potential to find accurate vegetation indices that correlates with visual BYD severity. Secondly, the visibility of BYD symptoms on canopies is likely to require finer pixel resolution to be reflected in digital images, whereas the imaging sensor size, the sensor total effective pixel resolution, and the UAS flight height adopted in this study might induce a ground sampling distance that could be further improved to discover clear BYD related features. Lastly, image acquisition in this study was based on a weekly frequency. The temporal resolution of data collection may not be sufficient to match the optimum period for observation of BYD symptoms. Disease scoring using HTP is scaling fast among breeding programs; however, how to effectively use this data remains challenging. Based on the variable heritability observed for BYD severity across seasons (0.211–0.806), genetic progress based only on phenotypic selection will be limited. In a year when there is a BYD outbreak, HTP can rapidly provide quantitative measurements compared to the alternative visual breeder score. Previous studies have shown that data collected with sensor-based HTP can be substituted to improve conventional disease visual evaluation (Sankaran et al. 2010; Kumar et al. 2016; Zheng et al. 2018); although our study is the first attempt to characterize BYD severity in wheat using HTP.

Genome-wide association analysis

Using GWAS we detected QTLs on chromosomes 5AS, 7AL, and 7DL for BYD severity BLUPs values. Using GBS tags that mapped to known alien fragments, we confirmed that Bdv2 resistance gene was located at 7DL and confirmed that the 7DL QTL was explained by the presence of the Bdv2 resistance gene. Even though only 33 wheat lines were positive for the presence of Bdv2, we still had enough power to detect its effect, supporting that Bdv2 has a strong effect on BYD under Kansas field conditions (Fig. 5). The associations on chromosome 7AL, observed for both BYD severity and Bdv2, suggest that the SNP markers on the 7AL peak may be miss-anchored markers that should have mapped to 7DL. The relatively high heritability values obtained for the untreated replications (Fig. 2) allowed us to detect a minor QTL on 5AS. Marza et al. (2006) reported a QTL at 38cM on the short arm of chromosome 5A associated with yellowing symptoms caused by BYD, and it is possible that this is the same region yet more data is needed to confirm if these QTLs are the same. The only other study reporting GWAS for BYD in wheat was able to identify several markers associated with BYD resistance on chromosomes 2A, 2B, 6A, and 7A (Choudhury et al. 2019b). However, most of the association were explained by individual SNP markers, and to date do not have any definitive biological link. GWAS results for the other traits used in this study did not discover genomic regions associated with the traits (Supplementary Fig. 4). Taken together, these results suggest that BYD resistance in the primary pool of wheat is rare and there is limited large effect loci that could easily be incorporated into the breeding program, thus GS could be an efficient way to enhance BYD resistance.

Genomic selection

We evaluated several different GS models to identify the best approach for predicting BYD (Fig. 6). Overall, we observed some trends including (1) incorporating years with consistent BYD disease data in the training population increased the model predictive ability; (2) predicting years with high disease pressure is difficult; and (3) using Bdv2 as a covariate had increased prediction performance, suggesting that it is responsible for much of the predictive power. These results suggest that GS based on G‐BLUP with Bdv2 as fixed effect covariate would lead to the greatest genetic gain for BYD breeding. Using selected major QTL as a fixed effect to improve GS models was suggested in a simulation study (Bernardo 2014) and demonstrated with empirical studies (Rutkoski et al. 2014). Nonetheless, using Bdv2 as a fixed effect covariate in our GS strategies did not consistently improve the predictive ability for PTHT_M or reduction in GY (Rice and Lipka 2019). However, there was not a consistent distribution of Bdv2 allele across the cohorts. GS predictive abilities for BYD were low compared to other disease (reviewed by Poland and Rutkoski 2016). However, since this is the first report of GS for BYD resistance in wheat, we do not have similar results to make better comparisons. One possible explanation we did not explore is if these lines were selected earlier in the breeding pipeline for BYD and therefore represent a poor training population for testing GS models. Another way to improve the predictive ability could be using multitrait GS models. There are some examples in the literature where using correlated traits to the trait of interest resulted in higher genomic prediction accuracies (Jia and Jannink 2012; Rutkoski et al. 2016; Crain et al. 2018). BYD has traditionally been reported to have low $H^2$ (Tola and Kronstad 1984; Choudhury et al. 2019b) and in this study, even with well managed plots that often had $H^2$ approaching 0.8, we still had difficulty reproducing these results year to year as evidence of the challenge of studying this pathosystem. Moreover, the correlation between HTP parameters and BYD phenotypes was interesting, but not sufficient to be useful in combination with GS in the germplasm tested.

Conclusions

We were able to show that Bdv2 has a major effect controlling BYD resistance in the KSU breeding germplasm. Apart from the known Bdv2 and a potentially novel 5AS region, we did not find evidence of other regions controlling BYD resistance supporting the hypothesis of limited resistance available in the current wheat gene pool and the highly polygenic nature of the trait. Moreover, our study was the first attempt to characterize and improve BYD field-phenotyping using HTP and apply GS to predict the disease. HTP traits showed strong correlation patterns with BYD severity, however, none of these parameters shared a common genetic architecture with BYD severity. The GS predictive ability results that we found in this study open the door for further improvement and testing GS implementation for breeding for BYD resistance. Continuing the improvement of BYD characterization and the search of new sources of resistance using species related to wheat, will be crucial to broadening the resistant genes available to introgress into wheat germplasm.

Data availability

Supplementary material, including raw and analyzed phenotypic data, genotypic data, and basic plot scripts are available at Dyrad doi:10.5061/dryad.ncjsxkswd and GitHub https://github.com/umngao/wsm1_bdv2.

Supplemental material is available at G3 online.

Funding

This material is based upon work supported by Kansas Wheat Commission Award No.: B65336 “Integrative and Innovative Approaches to Diminish Barley Yellow Dwarf Epidemics Kansas Wheat.” PS was supported through a US Fulbright-ANII Uruguay Scholarship. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of industry partners.

Conflicts of interest

The authors declare that they have no conflict of interest.