Screening of stable resistant accessions and identification of resistance loci to Barley yellow mosaic virus disease

Plant material

A collection of 181 barley accessions sourced from around the world were used for GWAS in this study (Table S1). All the accessions were firstly genotyped using the genotyping by sequencing (GBS) approach. Over 17,000 SNP markers were scored with the assistance of USDA-ARS/Kansas State University (Fan et al., 2017). After omitting those markers with a missing rate ≥20%, a minor allele frequency (MAF) ≤0.05 and heterozygosity ≥5%, 3,818 SNP markers with unique physical position (Based on reference genome Morex V3, Mascher et al., 2021) were obtained. Markers associated with known BaYMV disease resistance genes were collected from previous studies (Table S2), and these markers were further genotyped in these 181 accessions using a standard PCR approach (Shi et al., 2019). Afterwards, 21 gene markers with polymorphism were confirmed and added to the total molecular markers (Table S2). Therefore, a total number of 3,839 molecular markers (Tables S3 and S4) were used for linkage disequilibrium (LD) analysis and GWAS.

Infected leaves of two varieties, Supi 1 (only susceptible to BaYMV) and Dan 2 (only susceptible to BaMMV), were collected to identify the differences in sequences between BaYMV and BaMMV.

Field experiments

All plant materials reported in this study were autumn-sown (late October) and cultivated at the BaYMV disease identification nursery (BaYMV/BaMMV-infected) and via the use of three biological replicates. Two nursery sites were selected in this study. One site is located in the main campus of Yangzhou University (32°23′36″ N, 119°25′36″ E) with an average rainfall of 63.76 mm and daily temperature ranging from 6 °C to 15 °C and the other one is located in Yancheng (Jiangsu Institute for Seaside Agricultural Sciences, 33°25′17″ N, 120°12′ 48″ E) with an averaged rainfall of 40.52 mm and daily temperature ranging from 4 °C to 15 °C. Disease symptoms were scored in March. The experiments in Yangzhou were conducted in 2015, 2016, and 2020 growing seasons (2015 Yangzhou, 2016 Yangzhou, and 2020 Yangzhou), and those in Yancheng were in 2019 and 2020 growing seasons (2019 Yancheng and 2020 Yancheng). All accessions were cultivated under a conventional water and fertilizer management regime.

Visual evaluation of BaYMV disease resistance

For obtaining accurate visual evaluation, disease severity scores were conducted four times within seven-day intervals. Severity scores of 0 to 4 represent the damage caused by virus infection from light to severe (Fig. 1): 0 means healthy leaves with no symptoms of disease, 1 represents typical chlorotic spots, 2 indicates a quarter of the leaves display a yellow mosaic pattern, 3 represents more than half of the leaves display a yellow mosaic pattern, and 4 indicates more than three-quarters of the leaves display yellow mosaic pattern (Pan et al., 2021). Afterwards, a standardized area under the disease progress steps (sAUDPS) index was used for the evaluation of disease resistance, which can effectively combine and transform multiple investigations of disease grades into a single value for further GWAS analysis (Simko & Piepho, 2012).

Data analysis

Descriptive statistics, analysis of variance (ANOVA), Pearson correlation analysis, and t-test were analyzed using IBM SPSS Statistics 22 software. The lme4 package of R (Bates et al., 2015) was used to perform the best linear unbiased prediction (BLUP) of the sAUDPS index with repeats for different years (Fernando & Grossman, 1989). BLUPs of all the accessions were used for GWAS.

To analyze the stability of phenotypic values in different environments, the sAUDPS indices of the natural population at the Yangzhou site for 3 years and the Yancheng site for 2 years were analyzed by additive main effects and multiplicative interaction (AMMI) analysis. Firstly, principal component analysis (PCA) was used to analyze the interaction between the sAUDPS index and the environment, and the significance of each principal component was judged. If the principal components were significantly different, it means that they meet the requirements of model analysis (Crossa et al., 1991; Lei & Xi, 2006; Jiwuba et al., 2020). The AMMI analysis was analyzed by using the Agricola package in R (https://rdrr.io/cran/agricolae/man/agricolae-package.html). The equations of the model and the calculation formulas of the stability parameters are as follows:

$$y_{ijk}=\mathrm{\mu }+{\mathrm{\alpha }}_i+{\mathrm{\beta }}_j+\sum _{r=1}^n{\mathrm{\theta }}_r{\mathrm{\gamma }}_{ir}{\mathrm{\delta }}_{jr}+{\mathrm{\rho }}_{ij}+{\mathrm{\epsilon }}_{ijk}$$

$$D=\sqrt{\sum _{r=1}^N{W}_r{\mathrm{\gamma }}_{\mathrm{i}r}^2}$$

In these two formulas, $y_{ijk}$ is the trait of the i– accession in the j— environment the phenotypic value repeated for the K μtime, μ is the average value of the phenotype, ${\mathrm{\alpha }}_i$ is the genotype main effect of the i– accession, ${\mathrm{\beta }}_j$ is the environment of the j— environmental main effect, ${\mathrm{\theta }}_r$ is the r principal component eigenvalue under the interaction of genotype and environment, ${\mathrm{\gamma }}_{ir}$ is the genotype score of the r principal component under genotype by environment interaction, ${\mathrm{\delta }}_{jr}$ is the environmental score of the r principal component under genotype by environment interaction, ${\mathrm{\rho }}_{ij}$ is the residual error, and ${\mathrm{\epsilon }}_{ijk}$ is the test error. The stability parameter D value refers to the Euclidean distance between the position of an accession in the principal component space of the interaction and the origin. $W_r$ is the variance explanation rate of each principal component. The lower the D value, the more stable the phenotypic value.

The CMplot package of R (https://github.com/YinLiLin/CMplot) was used to plot our 3,839 markers and the distribution of known BaYMV disease resistance genes on chromosomes (Fig. 2). The LD analysis was analyzed as the factor that influences the power of detecting a QTL (Brito et al., 2011; Mucha, Bunger & Conington, 2015). The LD between pairs of SNPs was calculated as the squared allele frequency correlation (r²) using TASSEL 5.2 software (Bradbury et al., 2007) as described previously (Megerssa et al., 2021). The LD plot was drawn by OriginPro 2021b software (https://www.originlab.com/). In this study, when r² = 0.2, the LD of the association panel was approximately 894 kb (Fig. 3).

GWAS for BaYMV disease resistance

GWAS was performed using the rMVP package in R (Yin et al., 2021). Based on the observation of p-values in different models under quantile–quantile (Q–Q) plots, the FarmCPU was selected as the most suitable model (Kaler et al., 2017; Zhang et al., 2021). The Bonferroni of false discovery rate (FDR) correction of 5% calculated by the rMVP package was used for multiple comparison adjustment and as a threshold to identify significant marker-trait associations (Aguado et al., 2020; Megerssa et al., 2021).

In this study, the Bonferroni of FDR correction of the natural population was approximately 4.88. Markers at p ≤ 10^−4.88 (−log10(p) ≥ 4.88) level detected by GWAS were defined as significant marker-trait associations (MTAs).

The phenotypic variance explained (PVE) by MTA was calculated using the following formula:

$$\mathrm{P}\mathrm{V}\mathrm{E}={\displaystyle \frac{V_Q}{V_P}\times 100\mathrm{\%}}$$

$$V_Q=4f_{REF}{f}_{ALT}{\left[{\displaystyle \frac{1}{2}\left({\mu }_{REF}-{\mu }_{ALT}\right)}\right]}^2$$

In the formula, $V_Q$ represents the genotypic variance, $V_P$ represents the phenotypic variance, $f_{REF}$ and $f_{ALT}$ represent the frequency of the two genotypes of markers, respectively, and ${\mu }_{REF}$ and ${\mu }_{ALT}$ represent the average phenotypic values corresponding to the two genotypes (Wang, 2014).

The genes between the LD upstream and the LD downstream of the MTAs were therefore functionally annotated based on barley Morex V3 reference genome annotation (https://apex.ipk-gatersleben.de/apex/f?p=284:57).

Sequence analysis of BaYMV

At the height of the BaYMV disease period in Yangzhou and Yancheng, the infected leaves of susceptible barley accessions ‘Supi 1’ and ‘Dan 2’ were collected for RNA extraction. Total RNA was extracted by RNAiso plus (Takara Bio, Beijing, China) before the general quality and concentration were determined via the use of an OD-1000+ spectrophotometer™ (One Drop, Shanghai, China). Total cDNA was obtained by reverse transcription using a 1st strand cDNA synthesis kit (Takara Bio, Beijing, China). A total of 50–60 ng of 1st strand cDNA was input as the template to a 20 μL multi-RT-PCR reaction for identification of BaYMV in leaves displaying disease symptoms (Zou, Wang & Xu, 2016; Shi et al., 2019). The 277 bp and 648 bp fragments were detected in ‘Supi 1’, corresponding to the coding sequence of HvACTIN and BaYMV specific fragments, respectively. The 277 bp and 940 bp fragments were detected by multi-RT-PCR in ‘Dan 2’, corresponding to the coding sequence of HvACTIN and BaMMV-specific fragments, respectively (Fig. 4A). The results showed that ‘Supi 1’ was only infected with BaYMV and ‘Dan 2’ was only infected with BaMMV at the two assessed sites.

BaYMV and BaMMV terminal binding protein VPg encoding sequences were amplified by specific primers (Table S5) following a procedure similar to previously reported (Jiang et al., 2022). Targeted fragments were verified by Sanger sequencing (Beijing Genomics Institute, Beijing, China). The sequencing data were assembled using MEGA X software (Kumar et al., 2018) and saved as a FASTA file (Table S6).

Phenotypic variation of BaYMV disease

The average incidence of BaYMV disease indicated that the mildest outbreak occurred in Yangzhou in 2020 while the most severe outbreak occurred in Yancheng in 2019 (Fig. 5A). The sAUDPS indexes were significantly positively correlated (p < 0.01) in all environments (Table 1).

In the analysis of variance, accessions and sites were considered fixed and significant differences (p < 0.01) were found between accessions and sites for all dependent variables (Table 2). These results indicated that the incidence of BaYMV disease was not only affected by genetic factors, but also by environmental factors, the year of analysis and even ecological conditions from the testing sites.

Stability analysis of resistance to BaYMV disease

PCA was used to perform the interaction analysis between phenotype and environment. The result showed that there were significant differences (p < 0.01) among the values of Interaction Principal Components Axes 1, 2, and 3 (IPCA1, IPCA2, and IPCA3) (Table S1; Table 3), indicating that there is a significant interaction between the sAUDPS index and environment. The accumulative contribution rate of the three IPCA was 92.60%.

According to the significant principal component score, the stability parameter D value of the sAUDPS index of each accession was calculated and sorted in ascending order. The D value was expressed on the AMMI biplot by the distance between the accession to the origin (Fig. 5B). Among the top ten accessions, seven showed stable resistance to BaYMV disease, all having a lower sAUDPS index in different sites and years (Table 4).

Marker-trait associations

In this study, the FarmCPU approach was employed to detect the associations between the markers and BaYMV disease resistance. QLB1 (p = 1.29E−06) and SNP2571 (p = 6.94E−13) associated with the disease in Yangzhou were detected on chromosomes 3H at 612 Mbp and 6H at 557 Mbp, respectively (Fig. 6A). Among these MTAs, QLB1 explained 28.58% variation and SNP2571 explained 24.21% variation at the Yangzhou site, indicating a high cluster of mapped genes and a novel resistance locus, respectively (Table 5).

A total of five MTAs were detected in Yancheng. Three highly significant ones, QLB1 (p = 5.03E−06), SNP4036 (p = 1.81E−07), and SNP0024 (p = 6.66E−06) had PVE values ranging from 19.23% to 31.80% and they are located on chromosomes 3H at 612 Mbp, 5H at 40 Mbp and 7H at 18 Mbp, respectively. Two less significant MTAs, SNP3426 (p = 3.36E−06) and SNP1839 (p = 6.21E−06), had PVE values ranging from only 0.28% to 2.74% (Fig. 6C, Table 5).

Sequence analysis of VPg

The BaYMV VPg coding sequence (880 bp fragment) was amplified using ‘Supi 1’ collected from both sites, and the BaMMV VPg coding sequence (830 bp fragment) was amplified from ‘Dan 2’ collected from both sites (Fig. 4B). After sequencing, multi-sequence alignment was performed using the MUSCLE program in MEGA X software. There were twenty-two SNP differences in the coding sequence of YZ-BaYMV and YC-BaYMV, among which thirteen non-synonymous mutations resulted in eleven amino acid changes (Fig. 4C). The rate of amino acid mutation between two VPg coding sequences was 5.34% (11/206). There was no SNP difference between YZ-BaMMV and YC-BaMMV.

Detection of known resistance genes to BaYMV disease

In this study, markers associated with known BaYMV disease resistance genes were used for genotyping the natural population, and twelve markers showed significant associations (p < 0.01) with disease resistance (Table S2). Some known disease resistance genes derived from a single accession are rarely used in disease resistance breeding (Le Gouis et al., 2004; Kai et al., 2012; Johnston et al., 2015; Pidon et al., 2020). The MTA detected in multi-environments was QLB1, which was the diagnostic SSR marker for rym4/rym5 (Tyrka et al., 2008). In this study, the QLB1 showed the closest association with BaYMV disease resistance (Table S2), indicated by the highest PVE (31.80%) (Table 5). Therefore, the known BaYMV disease resistance genes rym4/rym5 are still the major resistance genes in the natural population used in this study.

Discovery of novel resistance loci to BaYMV disease

In this study, the FarmCPU was used to detect the associations between the molecular markers and BaYMV disease resistance under different BaYMV/BaMMV-infected nurseries. In Yangzhou, an MTA (SNP2571) which explained 24.21% of phenotypic variation was identified at the terminal region of chromosome 6H (Figure 6A), different from the known BaYMV disease resistance genes Rym14^Hb (Pidon et al., 2020) and rym15 (Le Gouis et al., 2004) which are mapped at the other terminal region of chromosome 6H (Fig. 2). Based on the LD decay distance, the genes between the 894 kb upstream and 894 kb downstream of SNP2571 were defined as candidate genes. According to the annotation information of the Morex V3 reference genome, a total of 14 typical plant disease resistance genes (R genes) were enriched in this region (Table S7). Among them, two genes encoding leucine-rich repeat proteins (HORVU.MOREX.r3.6HG0631890 and HORVU.MOREX.r3.6HG0632320), five genes encoding NBS-LRR disease resistance protein-like proteins (HORVU.MOREX.r3.6HG0631950, HORVU.MOREX.r3.6HG0632390, HORVU.MOREX.r3.6HG0632400, HORVU.MOREX.r3.6HG0632430 and HORVU.MOREX.r3.6HG0632630) and seven genes encoding homologs to the disease resistance protein RPM1 (HORVU.MOREX.r3.6HG0632440, HORVU.MOREX.r3.6HG0632480, HORVU.MOREX.r3.6HG0632490, HORVU.MOREX.r3.6HG0632530, HORVU.MOREX.r3.6HG0632550, HORVU.MOREX.r3.6HG0632610 and HORVU.MOREX.r3.6HG0632650) putatively represent the candidates genes for the novel BaMYV disease resistance QTL on chromosome 6H.

At present, only one BaYMV disease resistance gene, rym3/Undesigned has been identified on chromosome 5H in barley, which was mapped within a 174 Mbp range at a position 362–536 Mbp on the chromosome (Saeki et al., 1999; Werner et al., 2003b). In this study, one new MTA (SNP4036) was detected at the physical position of 40 Mbp on chromosome 5H (Table 5), which could explain 19.23% of phenotypic variation in Yancheng. Based on the LD decay (894 kb) of the association panel, a total of 23 annotated barley genes were identified for the disease resistance QTL (Table S7).

rym2 and rym7t mapped at the interval position from 213 Mbp to 369 Mbp on chromosome 7H showed resistance to BaYMV disease (Gotz & Friedt, 1993; Takata et al., 2012). In this study, SNP0024 at a position 18 Mb from the start of chromosome 7H was significantly associated with BaYMV disease. This novel and significant MTA explained 19.79% of phenotypic variation in Yancheng (Table 5). There were 81 annotated barley genes in the confidence interval of the novel resistance QTL near SNP0024, of which two genes encoding NBS-LRR disease resistance like proteins (HORVU.MOREX.r3.7HG0644180 and HORVU.MOREX.r3.7HG0644290) and another two genes encoding homologs of the disease resistance protein RPM1 (HORVU.MOREX.r3.7HG0644780 and HORVU.MOREX.r3.7HG0644810) were potential candidate genes for the novel disease resistance QTL on chromosome 7H (Table S7).

Differences of BaYMV disease in different nurseries

There have been many reported disease resistance genes that maintain variable levels of resistance to diverse BaYMV viral strains (Zheng et al., 1999; Jiang et al., 2020). In this study, variance analysis showed that there were significant differences (p < 0.01) in the natural population of BaYMV disease between these two sites. The mean sAUDPS index in Yangzhou from 3 years was lower than that in Yancheng from 2 years, indicating that the difference observed was due to differences in viral pathogenicity. Moreover, each testing site had independent resistance loci, except for QLB1 which was detected at both sites. Previous studies have shown that there were differences between strains of BaYMV in Yangzhou and Yancheng (Chen et al., 1999). Our results showed eleven amino acid differences in the BaYMV-VPg coding sequence between Yangzhou and Yancheng. VPg is a multifunctional protein, which directly participates in almost every key process of virus infection via the formation of different splice forms. VPg especially functions in protein maturation and virus movement during the process of infection, thus fundamentally determines viral pathogenicity (Stein et al., 2005; Charron et al., 2008; Li et al., 2016). Therefore, the difference in the VPg coding sequence of the BaYMV virus might be the reason for the difference in the severity of the BaYMV disease observed across the two testing sites.

Prospects of BaYMV disease resistance breeding

Under the natural environment, the dormant spores of Polymyxa graminis carrying virus particles exist in the soil and diseased root residues. Cultivation measures such as applying sterilizing agents, crop rotation and postponing the sowing date may alleviate the effect of BaYMV disease but these measures are not ecologically nor environmentally friendly or sustainable. In this way, the most efficient and environment-friendly way to alleviate BaYMV infection is to cultivate varieties of resistant species.

The relationship between virus and host plant disease resistance follows a typical pattern of gene-gene interaction. Although there are many identified genes regulating crop resistance to BaYMV disease, however, due to the rapid evolution of the virus, the positive effects from these genes are becoming increasingly diminished. Identifying novel genes conferring enhanced crop resistance to the evolving BaYMV is thus very necessary (Yang et al., 2014a; Yang et al., 2017) for offering new genetic resources for breeding accessions of barley resistant to BaMYV disease (Yang et al., 2017; Jiang et al., 2020).

In this study, the MTA (SNP2571) detected on chromosome 6H consisted of two alleles (C-allele and T-allele). A t-test conducted on the phenotype between the two groups at the Yangzhou site showed that the disease symptoms of the accessions containing the C-allele were significantly (p < 0.01) lower than those of the accessions containing the T-allele. In our study, the disease resistance allele of SNP2571 was mostly found in six-rowed landraces (n = 59, 32.60%). Similarly, the resistance allele of the MTA on chromosome 7H was also mostly discovered in six-rowed landraces (n = 44, 24.31%) (Table S8), supporting previous reports that novel genetic resources of barley resistant to BaMYV disease derived from landraces (Shi et al., 2019; Jiang et al., 2020) could be used in plant breeding for disease resistance. It should be emphasized that many R genes were enriched in the candidate interval between SNP2571 and SNP0024 based on LD analysis (Table S7).