Received 6 December 2015; accepted 24 January 2016; published 27 January 2016
1. Introduction
Puccinia striiformis Westend. f. sp. tritici Erikss. (Pst) is an obligate biotrophic fungus causing stripe (yellow) rust mainly on wheat. Significant yield losses have been reported in most wheat production regions, and the pathogen still threatens wheat production worldwide [1] [2] . Intensive efforts have focused on monitoring the disease, identifying and incorporating resistance genes into wheat cultivars, and use of fungicides to manage the disease [3] [4] . However, no genes for pathogenicity or virulence have been characterized for Pst.
To understand the mechanism of the “boom and bust” cycle of the battle between avirulent/virulent pathogens and resistant/susceptible cultivars, in the middle 20th century Flor made a breakthrough by demonstrating the gene-for-gene concept that “for each R gene that conditions resistance in the host there is a corresponding Avr gene that conditions pathogenicity in the parasite” [5] . Since then, scientists have begun to search for the biochemical and molecular basis for the gene-for-gene concept, and the elicitor-receptor model has been proposed. Based on this model, the recognition of an elicitor (product of an Avr gene) in the pathogen by a receptor (product of the matching R gene) in the resistant host plant activates a signal transduction pathway that leads to a resistant response, such as a hypersensitive response (HR), in the host [6] . Interestingly, most of AVR proteins among filamentous pathogen effectors have a signal peptide region [7] . Particularly in flax rust, Melampsora lini, four putative Avr genes (AvrL567, AvrM, Avr4, and AvrP123) have been identified from 21 haustorially expressed secreted proteins (HESPs), and transient expression of these genes in flax causes resistance gene- mediated cell death, thus demonstrating their avirulence activities [8] . There are some exceptions. For example, two protein effectors in the stem rust fungus (P. graminis f. sp. tritici) may activate the R protein long before haustoria formation [9] . The Pst-wheat pathosystem well fits the “gene-for-gene” concept. Although numerous resistance genes in wheat have been identified with some cloned [10] , no avirulence genes have been molecularly characterized.
Limited studies were conducted to identify candidate avirulence gene (Avr) in Pst. For example, Huang et al. [11] studied the expression profiling of 442 gene transcripts of Pst in in vitro and Pst infected wheat leaves using a custom-made Pst gene chip including secreted protein (SP) genes. The differentially expressed transcripts were involved in various biological processes. Several transcripts were expressed differentially between compatible and incompatible interactions, and thus could be virulence/avirulence related genes. Using quantitative RT-PCR analysis, Cantu et al. [12] identified 22 effector candidate genes which were significantly expressed during the infection in wheat. These candidate genes encoded secreted proteins whose transcripts were enriched in haustoria. Similarly, potential Avr genes were identified in cDNA libraries made from germinating urediniospores and fungal cells harvested during the infection processes by examining expression profiles of putatively secreted protein genes [13] [14] . All of these studies to identify candidate Avr genes were based on gene expression levels in germinated urediniospores and in infected wheat tissues, while phenotypic virulence profiles were rarely considered. In the present study, we attempted to identify Avr genes, corresponding to resistance genes in wheat, through association analysis between genetic variation and virulence/avirulence phenotypes taking the advantages of SP-single nucleotide polymorphisms (SP-SNPs).
SP-SNPs derived from expressed sequence tags (ESTs) would be promising markers for identifying avirulence genes since some SP-SNP markers may be closely associated with or directly in the coding regions, and thus more likely to be involved in causal alleles or highly correlated with the causal alleles. It is possible that the effectors in Pst are also highly expressed secreted proteins (HESPs) as in the study with M. lini [8] , because M. lini and Pst are related rust fungal species, and both fungi are obligate biotrophic parasites and obtain their nutrient from living hosts through haustoria. Moreover, most of the identified filamentous pathogen effectors have a signal peptide domain [7] . This would be the first study to identify Avr genes in Pst using an association mapping approach based on SP-SNP markers, and may lead to a feasible approach for identification and characterization of Avr genes in Pst, which is genetically and molecularly hard to manipulate.
2. Materials and Methods
2.1. Isolate Selection
A total of 352 Pst isolates used in the present study were selected from the P. striiformis collections made in 2010 and 2011 in the US. These isolates were selected to cover as many races as possible, thus representing a diverse combination of avirulence/virulence profiles.
2.2. Virulence Tests
The Pst isolates were tested for their avirulence/virulence patterns using the 18 Yr single-gene lines and supplementary differentials, as described by Wan and Chen [15] . The avirulence or virulence of isolates to a particular Yr resistance gene line was represented by an infection type (IT), which was scored using a 0-to-9 scale, with 0 as the most avirulent and 9 as the most virulent, as described previously [16] . ITs were recorded at 20 to 22 days after inoculation. For association analysis in this study, a phenotypic trait of an isolate to a particular Yr gene was defined as avirulence when IT was 0 - 4 and virulent when IT was 5 - 9. Totally, 19 phenotypic traits (18 Yr single-gene differentials plus one supplementary differential supposed to have a single stripe rust resistance gene, Yr2) were analyzed for associations with SP-SNP markers. These 19 resistance genes were: Yr1, Yr2, Yr5, Yr6, Yr7, Yr8, Yr9, Yr10, Yr15, Yr17, Y24, Yr27, Yr32, Yr43, Yr44, YrExp2, YrSP, YrTr1, and YrTye. Their corresponding avirulence genes were symbolized as AvYr1, AvYr2, and so on.
2.3. DNA Extraction
Genomic DNA was extracted using a modified CTAB method, as described previously [17] . The DNA stock solutions were diluted to 5 ng/µL with sterilized ddH2O. The DNA quality and concentration of the diluted working solution were checked by PCR amplification using the internal transcribed spacer (ITS) primers ITS3 and ITS4. The quantified working solution was used for SP-SNP genotyping.
2.4. Isolates Genotyping
Genotyping of 352 Pst isolates over SP-SNP loci was conducted on a Sequenom MassARRAY iPlex platform. In the initial PCR reaction, the locus-specific region about 100 bp containing a SN-SNP site was amplified. Then shrimp alkaline phosphatase (SAP) reaction was conducted in which phosphatase was used to eliminate unincorporated dNTPs from the previous reaction. A single base extension reaction was conducted using modified dideoxynucleotide terminators and thermosequenase enzyme, which produced two extension products differing only in the last one base for each SNP locus. The SP-SNP alleles were identified based on the distinct masses of the extended primers using the matrix-assisted laser desorption ionization-time-of-flight (MALDI- TOF) mass spectrometry. The automatically generated genotype calls were checked manually and overwritten based on the intensity peaks of the alleles in the mass spectrum using Spectro TYPER (SEQUENOM, Inc. San Diego, California). Generally, the locus having an intensity peak of one allele which was 2 times higher than the other one was considered as a homozygous locus, otherwise as heterozygous. Also, the background noise was taken into consideration for overwriting the automatic calls.
2.5. Data Analyses
The SP-SNP data were examined to determine the minor allele frequency (MAF), number of heterozygotes, and number of missing alleles. Correlations between SP-SNPs and the 19 virulence/avirulence traits of the Pst isolates were analyzed, and a Spearman correlation matrix was generated in the program Genome Association and Prediction Integrated Tool (GAPIT) [18] . The IT distribution of isolates of the 19 virulence/avirulence traits was also determined.
The mixed linear model (MLM) including both fixed and random effects was incorporated in the program GAPIT and used for association analyses. By including individuals as random effects, the MLM model was able to incorporate relationships among isolates in the association analysis. The relationships among isolates were represented as a kinship matrix computed using the VanRaden method. Specifically, two MLM scenarios were performed in this study, a naive scenario and a basic scenario with principal component analysis (PCA) and kinship. The naive scenario was achieved by setting to group using commend “group. from = 0 and group. to = 0”. By doing this, the scenario was less powerful by ignoring the fixed effects such as population structure, but served as a control for comparison with the basic scenario with PCA and kinship. While the basic scenario incorporated PCs as fixed effects and used to group kinship matrix for isolates relationships [18] . The commands were: “myGAPIT <- GAPIT(Y = myY, G = myG, group. to = 0, group. from = 0)” and “myGAPIT <- GAPIT(Y = myY, G = myG, PCA. total = 3, Model. selection = T, group. from = 351, group. to = 351)” for naive scenario and basic scenario with PCA and kinship, respectively. The false discovery rate-adjusted P values were applied to reduce the false-positive rate in the multiple testing carried out in the association analyses [19] .
3. Results
3.1. Infection Type Distribution
IT was used to qualitatively describe the virulence of Pst isolates. However, the IT data was not completely qualitative or quantitative, and no normal distribution was expected, as shown in Table 1. The balanced distribution of avirulence (e.g. IT 0 - 4) and virulence (e.g. IT 5 - 9) was expected for association and would avoid spurious association to some extent. For example, 193 out of 352 Pst isolates in the present study were highly avirulent (IT 1 - 2) to Yr1, and 157 were highly virulent to Yr1 (IT 8) while only 2 had an intermediate IT. In this case, the avirulence/virulence distribution for Yr1 was relatively balanced, nearly 1:1 for avirulent to virulent isolates. A similar case was observed for YrTye (Table 1). However for some virulence phenotypes, the avirulence/virulence distributions were not balanced. For example, all isolates used in this study and those collected thus far from the US were avirulent to Yr5 and Yr15 [1] [15] . Therefore, these two traits were not included for further analyses. Even though the virulence/avirulence distributions of 15 other traits were not exactly balanced, they were included in the association analysis as either avirulent or virulent isolates were more than 5%.
3.2. Allele Tests
The frequency distributions of minor alleles (MAF, minor allele frequency), heterozygous, and missing alleles were shown in Figure 1. The mean frequency of minor alleles was around 40%. The Pst population used in the present study had a high heterozygosity (data not shown), but largely within an acceptable range. All alleles with a frequency less than 0.05 were filtered prior to analysis. There were very few missing data points. Generally, these estimates indicated that the genotypes of the Pst population in this study were suitable for further analyses.
3.3. Principal Component Analysis (PCA) and Kinship Analysis
Principal component analysis conducted using the GAPIT program indicated that the first two PCs explained the largest portions of the genetic variance, as shown in Figure 2. The Efficient Mixed Model Association (EMMA) algorithm was used to establish a kinship matrix, a heat map of values in the kinship matrix is shown in Figure 3.
![]()
Table 1. Distribution of infection types of 352 Puccinia striiformis f. sp. tritici isolates scored on 19 wheat lines with Yr single-genes.
![]()
Figure 1. Boxplots of the frequencies of minor alleles and heterozygous loci and missing data for 352 Puccinia striiformis f. sp. tritici isolates at 97 SP-SNP loci.
![]()
![]()
(a) (b)
Figure 2. Principal component analysis (PCA) of 352 Puccinia striiformis f. sp. tritici isolates using 97 SP-SNP markers. (a) The genetic variance explained by each principal component, up to 10. (b) Plot of the second principal component (PC2) against the first principal component (PC1).
Generally, two lineages were estimated for the Pst population in this study. Within each lineage, the isolates had a closer kinship to other isolates, indicating a strong population structure or stratification.
3.4. Thirty SP-SNPs Significantly Associated with Nine Pst Virulence Phenotypes
The two scenarios for association analyses in this study were compared in terms of quantile-quantile (QQ)-plot of P values and the optimum compressions. AvYr1 was taken as an example to illustrate the comparison between the two scenarios. In the QQ plot, the observed −log10(P) values were plotted against the expected −log10(P) values. Thus, the SP-SNPs under the null hypothesis of no association between the SP-SNP and the trait would be laid on the diagonal line in the QQ-plot. The SP-SNPs deviated from the diagonal line indicated potential association, either spurious due to population structure and isolates relatedness, or truly associated with the trait. The comparison of two scenarios with AvYr1 as an example was shown in Figure 4. Compared to naive scenario which did not take PCA and kinship into account for association analysis, much fewer SP-SNPs were deviated from the null hypothesis of no association in the basic scenario, indicating the spurious associations in the naive scenario. A low level of inflation of the test statistics in the basic scenario with PCA and kinship indicated a high efficiency in controlling spurious associations due potentially to the population structure and cryptic relatedness among the isolates. Similarly, for comparison of the optimum compression for AvYr1 using these two scenarios, a low (−2 × log) likelihood function in the basic scenario with PCA and kinship (Figure 5) with all Pst isolates grouped into 344 groups for kinship analysis, suggesting a high likelihood function, while only 4 groups
![]()
Figure 3. Kinship plot of 352 Puccinia striiformis f. sp. tritici isolates. The heat map was generated based on the values in the kinship matrix which was computed using the efficient mixed model association (EMMA) algorithm.
![]()
Figure 4. Comparison of quantile-quantile (QQ)-plot of P-values using different association models for the AvYr1 locus. In the QQ-plot, the observed negative base 10 logarithm of the P-values were against the expected observed negative base 10 logarithm of the P-values under assumption that the P-values follow a uniform [0,1] distribution. The dotted lines represented the 95% confidence interval for the QQ-plot under the null hypothesis of no association between the SP-SNP and the AvYr locus. (a) Computed using the naive scenario in the GAPIT program; (b) Computed using the basic scenario with PCA and kinship.
![]()
Figure 5. Comparison of the profile for the optimum compression using different association models for AvYr1. (Left) Computed using the naive scenario method. The optimal method to estimate group kinship was “Mean”, the optimum clustering method was “average”, the number of groups was 4, the value of −2 × log likelihood function was 1454.34, and the genetic variance (heritability) was 0.721. (Right) Computed using the basic scenario method with PCA and kinship. The optimal method to estimate group kinship was “Mean”, the optimum clustering method was “average”, the number of groups was 344, the value of −2 × log likelihood function was 1406.60, and the genetic variance (heritability) was 0.402.
were accounted for kinship. Evidence that the basic scenario with PCA and kinship was superior to the naive scenario model was indicated by comparisons of these two scenarios for other virulence traits (data not shown).
The mixed linear model using the basic scenario with PCA and kinship indicated associations between 30 SP-SNP markers and 9 avirulence genes (corresponding to Yr genes) with a false discovery rate (FDR)-adjusted P value lower than 0.05 (Table 2). These 9 genes were AvYr1, AvYr2, AvYr7, AvYr8, AvYr44, AvYrExp2, AvYrSP, and AvYrTye. For all these avirulence genes, each was significantly associated with two or more SP-SNPs, except AvYr2 and AvYrExp2. Also, some SP-SNPs were significantly associated with two or more avirulence genes. For example, significant associations were indicated between SP-SNP marker PstSNP00438 and AvYr1, AvYr7, and AvYr44. Specifically for AvYrTye, all three significantly associated SP-SNPs were also significantly associated with AvYr1. This result was not surprising since the virulence ITs of these traits were not completely independent. To confirm this, the correlation relationships among the 19 traits were estimated, and the results were given in Figure 6. The correlations between most of the avirulence genes were low or moderate in terms of correlation coefficients. However, high correlations were suggested between AvYr1 and AvYrTye (R = 0.90), and between AvYr7 and AvYr44 (R = 0.89).
4. Discussion
The preliminary results in the present study showed that association analysis can be a powerful tool for relating molecular markers to avirulence genes of Pst. Compared with previous gene expression studies [11] -[14] , the association analysis takes the virulence phenotypes into consideration, instead of focusing only on differentiation of gene expression of the pathogen in the Pst-wheat interaction process. Therefore, it may provide direct genetic evidence for virulence variation. Even though only 97 SP-SNPs were used in this study, we found significant associations (P_FDR_adj < 0.05) for 41 SP-SNP Pst virulence associations involved in 9 AvYr genes and 30 SP-SNP loci (Table 2). Most of the genes used in this study for SP-SNP development were studied previously by Huang et al. [11] who analyzed the gene expression profiling from both compatible and incompatible interactions in which “Avocet Susceptible” and resistant “AvSYr5NIL” were inoculated with race PST-78. Six of SP genes identified in the present study were also differentially expressed in their studies. These genes are
![]()
Table 2. SP-SNPs associated with avirulence genes of Puccinia striiformis f. sp. tritici.
a. Avirulence genes (AvYr) are corresponding to wheat Yr resistance genes. b. P_FDR_adj (false discovery rate-adjusted P values) were estimated from association analysis using basic scenario with PCA and kinship model from GAPIT program. c. Supercontig, the SNP locus position in the supercontig, and protein ID are from the whole genome sequence of PST-78 in BROAD Institute Puccinia database (http://www.broadinstitute.org/). d. Protein function is based on BLAST search in NCBI database. e. NS = non-synonymous; S = synonymous; and UTR = untranslated region.
![]()
Figure 6. Correlation coefficients between 19 avirulence/virulence loci of Puccinia striiformis f. sp. tritici. The avirulence genes are represented by their corresponding resistance Yr genes to save space.
represented by SP-SNPs 00021, 00024, 00363, 00392, 00414, and 00441. Interestingly, these genes were from either the urediniospore cDNA library or germinated urediniospore cDNA library, but none were from the haustorial cDNA library which was thought as an Avr gene reservoir in rust fungi [8] . Most of the identified genes in the present study encode hypothetical proteins with no known functions (Table 2), while others are involved in different biological processes, but none of them are known to be involved in fungal pathogenicity. One gene associated with AvYr1 and AvYrTye has a high homology with the mago nashi gene in Drosophila. This is an essential gene required for the assembly of germ plasm in Drosophila, and nonsense mutations in this gene caused zygotic lethality [20] . A glycerol kinase gene was identified as associated with AvYr6. Glycerol kinase has an activity of catalyzing the phosphorylation of glycerol to glycerol-3-phosphate (G3P), which is an energy-pro- ducing reaction. However, G3P contributes to basal resistance in plants against Colletotrichum higginsianum, a hemibiotrophic pathogen [21] . Another gene associated with AvYr6 encodes a Sm-like protein, LSm2 which binds to the U6 small nuclear ribonucleoprotein and may be involved in pre-mRNA splicing [22] . Also associated with AvYr6 is an alpha-trehalase gene which is involved in hydrolysis of trehalose, an important storage compound in vegetative cells and spores of fungi [23] . A gene encoding GDP-mannose 4,6-dehydratase was associated with AvYr7 and AvYr44. This enzyme is involved in bio-synthesis of GDP-L-fucose which is important to the development and strength of stem tissues in plants [24] . A high-mobility-group box1 (HMGB1)-coding gene was identified as associated with AvYr7 and AvYr44. HMGB1 has multiple functions in DNA architecture and transcriptional regulation. It activates immune cells involved in immune processes [25] . A gene associated with AvYrExp2 encodes a histone H4, which is a basic component of the nucleosome. Histone H4 is subject to covalent modification (e.g. acetylation, methylation), therefore it may change the expression of genes that are associated with its parent histone octamer [26] . We need answers for several questions before we can make clear conclusions about the associations of the SP genes and avirulence loci.
The first question is why not all heritability can be explained by the associated SNP markers. Even though significant associations were detected, all associated SP-SNPs accounted for only a limited amount of heritability ranging from 0.00 to 0.63 with an average of 0.26, which was consistent to most genome-wide association studies (GWAS) [27] -[29] . The major reason underlies the hypothesis that common diseases may be caused by common genetic variants. Therefore, in the association analysis the rare variants (for example MAF <5%) are not taken into consideration. However, these rare variants probably also cause common diseases by either acting alone or accumulating effects from many alleles [28] . MAFs of greater than 1% were proposed for SNPs in the GWAS [30] . Particularly in the present study, the SP-SNPs were filtered by MAF > 5%. It is reasonable to speculate that the rare variants may be responsible for the unbalanced virulence distributions which were common for most virulence loci (Table 1). Another possible reason for the unexplained heritability is that only SNP variation is evaluated in most association analyses including the present study. Even though SNP is a common form of variation in most organisms and is a causal variant for some diseases [31] , other variations (e.g. copy number variants) should also be incorporated in a GWAS [28] .
The second question is what genes they are and what functions they have. For the genes identified in this study, eight were annotated based on homology searches to known genes in the NCBI database. However, none of these genes, as far as we learn from the database, were directly involved in the pathogenicity in the fungi. Such non-causative variations identified in GWAS are not uncommon. A more complicated case is that a lot of markers identified from association analyses are not even in the gene regions. The non-causative associations, probably spurious associations, are mainly due to the linkage disequilibrium (LD) between many loci in the genome [30] [32] . In the LD, the non-causative variants and the causative variants are non-randomly associated, which commonly occurs in structured populations [33] . However, in this study we were unable to estimate the LD decay because of the limited number of SP-SNPs and unknown physical distances between them.
The third question is how to address false-positive associations in an association study, especially in which structured populations are used, like in this study. One source of false-positive associations is the statistical fluctuations arisen from chances from multiple-hypothesis testing [32] . To address this problem, some statistical approaches have been proposed including the Bonferroni correction and estimation of the FDR that are two common correction methods [34] . Another source of false-positive associations is from the structured populations in the association studies [32] . If the individuals from subgroups from such a population happen to differ in allele frequency of genetic polymorphisms which are unrelated to the trait of interest, the false-positive association will occur. The population structure, or population stratification, is expected in the Pst population because of the mainly clonal reproduction. Two major subgroups were identified, one avirulent to Yr1 and YrTye, with the other virulent to these genes. Therefore, any SP-SNPs that were different in allele frequency between these two subgroups might appear to be associated to AvYr1 and AvYrTye. To better rule out the spurious associations caused by population structures in Pst, sexually reproduced segregating populations should be developed, and the successful cases in plants, such as in maize, provided an example [29] [35] . To control the effects of population stratification for GWAS and take the advantage of linkage-association analyses, 25 diverse inbred maize lines were crossed to a reference line to capture a large number of recombination events, and 15 RIL families were produced with 200 RILs for each family [35] . This approach was called nested association mapping. Besides these two approaches, replicating genotype-phenotype associations in larger and independent populations was another way for establishing the credibility of GWAS [36] .
Fortunately, these obstacles in the GWAS to identify virulence or pathogenicity related genes in the Pst population can be conquered, to some extent at least, if not wholly. Firstly, a larger number of isolates should be selected and genotyped through whole genome sequencing approaches. With a large number of isolates, the rare variants are more likely to be detected in the population. The SP-SNPs at the genome coverage should be developed and used for GWAS. The SP-SNP approach will continue focusing on functional genes, instead of non-gene regions. High resolution mapping of an association study depends on the number of markers as well as on the LD decay. Therefore, the LD decay of the Pst genome also should be estimated. To develop more SP-SNPs in Pst, a project on deep sequencing 12 isolates with different virulence spectra is underway. These deep genomic sequences will enable us to develop more markers for GWAS, including gene absence/presence and deletion/insertion polymorphisms which are probably important in the evolution of virulence of Avr genes in plant pathogens. For example, novel avirulence genes in Magnaporthe oryzae, the rice blast fungus, were identified from absence regions in the genome [37] . Secondly, the recent discoveries of alternate hosts, on which sexual reproduction of Pst can occur [38] -[40] enable us to establish a mapping population. Such sexual Pst populations should provide an efficient way to control population stratification by breaking down the LD in the Pst genome. Also, with hundreds of thousands of genetic markers, estimation of a linkage map with a high resolution for avirulence genes in Pst is possible. The identification of Avr genes in the Pst will not be the ultimate goal, but the first step to study the function, the evolution, and the genetics of the Avr genes, and thus to study the molecular mechanisms underlying the interactions between Pst and its hosts.
Acknowledgements
This research was supported by the U.S. Department of Agriculture, Agricultural Research Service (Project No. 5348-22000-015-00D) and Washington State University (13C-3061-5665) PPNS No. 0707, Department of Plant Pathology, College of Agricultural, Human, and Natural Resource Sciences, Agricultural Research Center, HATCH Project Number WNP00663, Washington State University, Pullman, WA 99164-6430, USA. The Chinese Scholarship Council (CSC) scholarship given to Chongjing Xia is highly appreciated. We are grateful to Drs. Scot Hulbert and Tobin Peever for critical review of the manuscript.
NOTES
*Corresponding author.