In-depth secretome analysis of Puccinia striiformis f. sp. tritici in infected wheat uncovers effector functions

Puccinia striiformis f. sp. tritici (Pst) is an obligate biotrophic fungus, the causative agent of wheat stripe (yellow) rust disease, a disease that can drastically reduce the production of wheat worldwide. The use of fungicides is a suitable solution for combating Pst, but using resistant wheat strains to control Pst genetically would be better for the environment and more cost-effective. In order to achieve host resistance against rapid evolving virulent strains of Pst, understanding the molecular basis of the disease is crucial.

Pst infects the leaf of the wheat plant by penetrating to the interior through the stomata and later differentiating into various developmental entities before finally entering into the mesophyll cells in which it can generate a feeding structure called a haustorium. The main interface for pathogen–host interactions is in the haustorium that expands between the plant cell wall and the cell membrane. This interface can activate primary host defenses when pathogen-associated molecular patterns (PAMPs) are detected by pattern recognition receptors (PRRs), resulting in PAMP-triggered immunity (PTI) [1]. Virulence is achieved by Pst through the deployment of secreted proteins, termed as effectors, into the apoplastic fluid or inside the host [2]. The effectors are useful weapons in the fungal arsenal that allow Pst (i) to overcome PTI, (ii) to establish a suitable environment, and (iii) to help proliferation [3]. Particular effectors are recognized by cytoplasmic receptor proteins in the host cell, triggering a defense response called effector-triggered immunity (ETI). Compatible and incompatible interactions depend on whether the effectors can achieve virulence by evading Resistance (R) proteins, while the R proteins evolve to detect the effectors. Therefore, evolutionary pressure imparts a driving force on the race between host resistance factors and pathogen virulence factors [4].

Investigating the secreted effectorome of the fungus guides us toward a better comprehension of the nature of pathogenic virulence and determinants that can lead to avirulence. The combination of next-generation sequencing (NGS) and data mining approaches represent a great opportunity to determine candidate effectors secreted by the fungus. Several studies attempting to construct a secretome of either Pst or its closest relatives (stem rust pathogen Puccinia graminis f. sp. tritici (Pgt) and leaf rust pathogen Puccinia triticina (Ptt)) have been published. Yin and colleagues reported a set of secreted proteins generated from expressed sequence tag (EST) sequences of a cDNA library generated from Pst haustoria [5]. A repertoire of small-secreted proteins (SSPs) was reported using a pipeline for effector mining in the genome sequence information of Pgt and Melampsora larici populina (the pathogen that causes poplar leaf rust) [6]. The isolate PST-130 was sequenced and assembled to discover Pst genes and identify candidate effectors [7]. The genome sequences of Pgt (CRL 75-36-700-3), Ptt (BBDD), PST-78 (2K-041), PST-1 (3-5-79), PST-127 (08-220), PST-CYR-32 (09-001) were released by the Broad Institute (http://www.broadinstitute.org/). Cantu and colleagues used comprehensive genome analysis to identify candidate secreted effectors in two United Kingdom isolates (PST-87/7 and PST-08/21) and two United States isolates (PST-21 and PST-43) [8]. They also released data of infection-expressed genes (6 and 14 dpi) and haustoria-expressed genes (7 dpi) for PST-87/7. Another report was published using transcriptome sequencing on Pst isolate 104E137A to elucidate haustoria-secreted proteins (HSPs) and expose the best candidate effectors [9]. A pioneer correlation analysis between genomes of seven existing races and seven new races of Pst was pursued to predict avirulence (Avr) determinants among SSPs [10]. In this study, 7 new races of Pst were sequenced with NGS and then combined with 7 older published genomes for a total of 14 races of Pst that were then subjected to correlation analysis to point out Avr candidates [10].

Data mining methodology is a promising strategy to identify candidate-secreted effector proteins (CSEPs) in genome and transcriptome data, as the incorporation of newly generated data continually improves the prediction algorithms using deep machine learning. Likewise, our understanding of the infection and disease progression of the pathogen is steadily expanding. Therefore, new information (candidate effectors) can be generated by exhausting presently available sequence data. A combined effort of using both comparative transcriptomics and data mining to narrow down the ever-expanding sets of CSEPs will simplify and focus on the work being done to discover the function of each CSEP experimentally. Only a few reports have focused on integrating both in silico prediction and transcriptome sequencing data on particular Pst races at the defined time intervals of haustoria formation and infection [8,9].

In the present study, our objective is to identify differentially expressed gene sequences (DEGs) of Pst during the interaction with both susceptible and resistant host genotypes. YR10 gene is one of CC-NB-LRR type resistance (R) genes [11,12]. We acquired de novo transcriptome sequencing data from both susceptible (AvocetS) and resistant (AvocetYR10) wheat cultivars infected with Pst for 10 days in order to discover differentially expressed secretome repertoires composed of small secreted proteins at a unique time point of the infection process after haustoria formation. Detected PstDEGs were then subjected to in silico prediction analysis to identify differentially expressed small secreted proteins (PstDESSPs). Then, the cataloged secretome candidates were analyzed in-depth via annotation and prediction programs to sort them into classified functions. Furthermore, a comparison of published Pst expression data [8,9] indicated the presence of 14 overlapping mutual CSEPs as the most promising subset for experimental validation. To test the reliability of our work, cDNA of one of the CSEP genes (Unigene17495, homologous to PSTG10917) was generated from Pst-infected wheat leaves (10 dpi) and cloned into Agrobacterium tumefaciens-compatible expression vectors. Heterologous expression of the CSEP in Nicotiana benthamiana revealed it was localized to the chloroplast in plant cells. Furthermore, we observed the biological function of Unigene17495 to be suppression of the INF1-mediated cell death response to N. benthamiana.

Near-isogenic wheat lines of Avocet-S and Avocet-YR10 were used to study both compatible and incompatible interactions upon Pst inoculation. Pst isolates were collected from the fields of Anatolia, Turkey, and obtained from the ‘General Directorate of Agricultural Research and Policies (TAGEM) of The Ministry of Food, Agriculture, and Livestock’. Avocet-S cultivars are susceptible to PstTR (which has virulence toward Yr2, 6, 7, 8, 9, 25, A, EP, Vic, Mich), whereas Avocet-Yr10 is resistant. Seeds were planted in soil and grown in a growth chamber under ideal conditions (relative humidity 60%, 22°C, 16 h light, and 8 h dark). The Pst inoculations were performed when the seedlings reached the two-leaf stage by spraying freshly collected urediniospores directly on to the leaves. Before the inoculation, urediniospores (100 mg) were germinated for 5 min at 42°C and treated with mineral oil (1 ml) to enhance their attachment to the leaf blade. For the control samples (mock inoculations), only mineral oil was used. After air-drying for 20 min, the plants were kept in the dark at 10°C overnight in the presence of continuous mist in order for infection to take place. After overnight infection, the standard growth conditions were resumed. The wheat leaves were collected at 10 days post-inoculation (10 dpi) when the disease symptoms were visible on the susceptible cultivar.

The flow chart in Figure 1 describes the overall protocol followed for collecting infected leaf samples using homogenization with liquid nitrogen and a sterile mortar and pestle. For 0.1 g leaf sample, 0.02 g poly(vinylpolypyrrolidone) (PVPP) was added during mortar and pestle grinding. For each condition, three biological replicates were homogenized together. Total RNA isolation was carried out using QIAzol Lysis Reagent (Qiagen) by following the manufacturer’s protocol. A NanoDrop was used for the spectrophotometric quantification of the total RNA samples. The samples were shipped to Beijing Genomics Institute (BGI) in absolute ethanol for sequencing. BGI measured the RNA quality of the samples on an Agilent 2100 Bioanalyzer, and the samples with a RIN value above 8 were allowed to proceed to mRNA enrichment and cDNA synthesis using random hexamers. Single nucleotide adenine was added to the purified and end-repaired short cDNA fragments. After adapter ligation, size selection was used to obtain templates for PCR amplification. The Illumina HiSeq™ 2000 platform was used to perform de novo transcriptome sequencing. BGI-Shenzhen (Shenzhen, China) conducted all the protocol steps from mRNA enrichment to sequencing. All quality controls were evaluated during sequencing using the Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System for quantification and qualification of the sample libraries.

Before assembly, clean reads were obtained from raw reads using internal BGI software (filter_fq) that follows three main parameters: (i) elimination of reads with adaptors, (ii) removal of reads with unidentified nucleotides larger than 5%, (iii) elimination of low-quality reads (if the percentage of the low-quality value (Q ≤ 10) reads are higher than 20%, they are removed). For de novo transcriptome assembly, clean reads of each library were assembled into contigs and unigenes by the Trinity program (release-20130225, http://trinityrnaseq.sourceforge.net/) [13].

Unigenes were aligned using the Blastx program (e-value < 0.00001) and the following protein databases: NR (NCBI, non-redundant database) [14], Swiss-Prot (EMBL protein database) [15], KEGG (Kyoto Encyclopedia of Genes and Genomes) [16], and COG (Clusters of Orthologous Groups) [17]. The direction of the unigene sequences was determined using the best alignment scores. NR, Swiss-Prot, KEGG, and COG were sourced to predict the coding region (CDS) of unigenes in the subsequent order of priority if a conflict in sequence orientation arose. ESTScan software was performed on unigenes with no successive alignments [18]. The results for NR, NT, Swiss-Prot, KEGG, and COG database alignments were also used for unigene function annotation. To obtain gene ontology (GO) annotation from NR annotations, the Blast2GO program (v2.5.0) was performed (e-value < 0.00001) on fungi databases (taxa, 4751) [19]. Web Gene Ontology Annotation Plot (WEGO) software was used for the classification of GO annotations [20]. The KEGG database was employed to investigate the metabolic pathway analysis of unigenes [16].

We identified our DEGs through transcriptome data obtained from pathogen-inoculated wheat leaves. The project generated two types of data: one is a de novo assembled transcriptome of Pst. The other one is DEGs mapped to wheat D genome [23], which is advantageous to assess both compatible and incompatible interaction of wheat with Pst. We separated DEGs mapped to wheat (the data not shown), the results will be discussed elsewhere to analyzed host-associated changes against Pst. We proceeded with Avocet-Pst/Yr10-Pst combination to obtain de novo assembled transcriptome of Pst and focused on pathogen-mapped DEGs. However, we needed to eliminate any possible contaminating sequences from the host. To accomplish this, we constructed a local database (named ‘The Pucciniale database’) by merging protein data of Puccinia striiformis f. sp. tritici (isolate 2K41-Yr9, race PST-78), Puccinia graminis f. sp. tritici (strain CRL 75-36-700-3, race SCCL), and Puccinia triticina (isolate 1-1, race 1-bbbd). All protein data were obtained from the Broad Institute Puccinia website (http://www.broadinstitute.org/). We also merged protein data of the genome of Triticum aestivum (IWGCv.1, downloaded from ‘ensemble.plant.org’ website) [24] to ‘The Pucciniale database’. We performed a BlastX analysis (E-value cutoff ≤ 10⁻⁵ or 1e⁻⁵) on DEGs using our local database of Pucciniale and wheat [25]. If the hit matched with the wheat proteome or produced a significant hit with the Pucciniale proteome, we discarded the sequences following the decision rule. The best hits within the Pucciniale proteome were evaluated as PstDEGs. The sequences with no significant hits against neither the pathogen nor host were assessed as unknown novel DEGs.

We predicted the ORFs and protein sequences from the PstDEG data using the ORFPredictor web tool (http://bioinformatics.ysu.edu/tools/OrfPredictor.html). We followed a widely accepted secretome prediction pipeline to construct our own differentially expressed secretome database [6,26,27]. The predictions were made following three criteria: (i) presence of N-terminus signal peptide region, (ii) absence of transmembrane helices, and (iii) mature protein length should be smaller than 300 amino acids. SignalP (version 4.1) was used to predict the presence of a secretion signal [28]. We estimated transmembrane helices using the TMHMM webtool (version 2.0, http://www.cbs.dtu.dk/services/TMHMM/) [29]. Filtered sequences were termed as differentially expressed SSPs (DESSPs) of Puccinia striiformis f. sp. tritici (PstDESSPs). The same filtering pipeline was applied to total proteins obtained from the ‘Broad Institute database’ for Pst and Pst relatives (Pgt and Ptt) (http://www.broadinstitute.org/). Four secretome datasets (PstDESSPs, PstSSPs, PgtSSPs, PttSSPs) were subjected to effector prediction via the EffectorP 1.0 and 2.0 programs (http://effectorp.csiro.au/) [30,31]. The Blastp program was used (E-value < 10⁻⁵ or 1e⁻⁵) for characterization of PstDESSPs by scanning against the following databases: the lipase engineering database (http://www.led.uni-stuttgart.de/), the protease and protease inhibitors (MEROPs) database (http://merops.sanger.ac.uk/) [32], the carbohydrate-active enzymes (CAZymes) database (http://csbl.bmb.uga.edu/dbCAN/) [33], the fungal peroxidase database (fPoxDB, http://peroxidase.riceblast.snu.ac.kr/) [34], and the pathogen-host interaction database (http://phi-blast.phi-base.org/) [35]. The ClustalW program was executed for multiple alignment analysis of PstDESSPs (https://www.ebi.ac.uk/Tools/msa/clustalo/) [36]. The iTOL webtool was used to construct a phylogenetic tree using multiple alignment results (https://itol.embl.de/) [37]. Small secreted cysteine-rich effector proteins were assessed as cysteine-rich if the protein possessed four cysteine residues or more. A comparison between infection- and haustoria-expressed genes was made using a similarity search using Blastp (E-value cutoff ≤ 1e⁻⁵ or 10⁻⁵) on publicly available data [8,9]. We scanned for PstDESSPs overlapping with Avr determinant candidates originally interpreted using correlation analysis by Xia and colleagues [10]. Subcellular localization predictions were done on mature proteins (without N-terminus signal peptides) to forecast translocation sites in the host plant using four different tools: TargetP 1.1 (http://www.cbs.dtu.dk/services/TargetP/) [38], Localizer (http://localizer.csiro.au/) [39], WolfPSORT (https://wolfpsort.hgc.jp/) [40], and ApoplastP (http://apoplastp.csiro.au/) [41]. All the analyses, annotations, and predictions regarding PstDESSPs are displayed in Supplementary Table S3.

The cloning primers were designed using the Pstg10917 sequence with an N-terminus CACC site in forward primers (with and without signal peptide (SP)) and no stop codon in reverse primers (Supplementary Table S4). Total RNA samples of infected Avocet-S leaves (10 dpi) used in transcriptome sequencing were sourced to synthesize cDNA after DNase treatment. The Transcriptor First Strand cDNA Synthesis Kit (Roche) was used following the manufacturer’s protocol, except that we preferred using a mixture of both random hexamer primers and oligo(dT) primers at equal volumes to the reaction primer. The effector of interest, with and without SP, was amplified from cDNA, cloned into the pENTR/D-TOPO vector (Invitrogen), and then subsequently cloned into the pK7FWG2 expression vector [42] utilizing Gateway Cloning. Both colony PCR and sequencing analysis validated the cloning of constructs at the proper positions. Then, the plasmid constructs were electroporated into Agrobacterium tumefaciens GV3101.

A. tumefaciens infiltration experiments were performed using the previously described protocol [43] with slight alterations. A. tumefaciens GV3101 strains carrying expression vectors were plated on selective agar media. After 2 days, bacteria were scratched and suspended in sterile water. Bacteria were collected by centrifugation at 4000 rpm for 4 min at room temperature (24–25°C). Bacteria were washed two times with sterile water and two times with A. tumefaciens-mediated gene transfer-induction medium (10 mM MES, pH 5.6, 10 mM MgCl₂). The final cell concentration was adjusted to 0.4 at A_600nm for infiltration unless otherwise stated. The inoculums were injected into N. benthamiana leaves using a needleless syringe. The expressions of the gene of interests were observed in the next 2–4 days.

A. tumefaciens GV3101 cells carrying pK7FWG2/Pstg_10917 and pK7FWG2/ΔSP-Pstg_10917 were injected into N. benthamiana leaves (4–5 weeks old) for transient expression. After 2–3 days, leaves were cut into small pieces of 3–5 mm length and soaked in tap water. A Leica 385 TCS SP5 confocal microscope (Leica Microsystems, Germany) was used for subcellular localization imaging. The wavelengths used for imaging GFP fluorescence were excitation at 488 nm and emission at 495–550 nm. Alternatively, the excitation and emission wavelengths for chloroplast autofluorescence were in the far-infrared (>800 nm).

The pTRBO/GFP plasmid for GFP expression (courtesy of the Kamoun Lab, Sainsbury Laboratory, Norwich, U.K.) and the pK7FWG2/SP-GFP plasmid [44] for SP-GFP were used as negative controls in the cell death suppression assays. The A. tumefaciens GV3101 cells carrying the negative control plasmids (GFP and SP-GFP) or the effector of interest (with and without SP) were injected into a N. benthamiana leaf. After 24 h, A. tumefaciens GV3101 carrying cell death elicitor pGR106/Inf1 (a gift from the Bozkurt lab, Imperial College, London, U.K.) was injected in the same spot without damaging leaf tissue. Cell death was measured at 4 dpi. The infiltration ring area was assessed if the brown region was less than 25% considered as strong suppression, 25–50% as moderate, 50–75% as weak, and more than 75% as no suppression at all.

We have collected total RNA from Pst-infected wheat leaves of susceptible and resistant seedlings to analyze gene expression differences between the compatible and incompatible pathogen–host interactions. The comparative transcriptome strategy was employed in similar studies; for instance, Sonah and colleagues identified virulence factors of Leptosphaeria maculans on susceptible and resistant canola cultivars, whereas Gupta and colleagues, on the pea powdery mildew Erysiphe pisi pathogen of Medicago truncatula [45,46]. Similarly, we managed to gather data from Pst pathogen generated on both susceptible (AvocetS) and resistant (AvocetY10) genotypes. The sequence reads have been deposited in the NCBI BioProject database under the accession number PRJNA545454. A total of 9204786 clean reads for Pst-infected susceptible wheat (Avocet-S_PST) and 6899596 clean reads for Pst-infected resistant wheat cultivar (Avocet-YR10_PST) were obtained after quality filtering. Their Q20 values were 98.16 and 97.90%, respectively (Table 1).

The clean reads were assembled into more than 100000 contigs with an average length of greater than 200 nucleotides using Trinity software (Table 2). Unigenes were clustered based on their level of similarity with each other. A cluster was defined as containing unigenes with more than 70% similarity. We generated 12956 distinct unigene clusters for Avocet-YR10_PST and 13564 for Avocet-S_PST (Table 2). A total of 42275 and 55004 unigenes came from single genes within the two sets, respectively. These data are represented in Table 2 as the distinct singleton. The length distribution of assembled contigs and unigenes is provided in Supplementary Figure S1.

The unigenes obtained from both Avocet-S_PST and Avocet-YR10_PST were merged into a single integrated dataset to conduct an overall interpretation of transcriptome sequencing. The unigenes were aligned against the following databases: NR, NT, Swiss-Prot, KEGG, COG, and GO for annotation via the Blastx program (e-value < 0.00001). Overall, 26843 of 61105 unigenes, corresponding to 43.9%, were annotated in one or more databases (Table 3). The detailed annotation results are listed in Supplementary Table S1. The NR database provided the highest number of annotations. Approximately half of the annotated unigenes shared high similarity (>60%) and aligned to either Pst genes (41.4%) or genes of other fungal species (Figure 2). The results of COG classification and GO functional annotation of the unigenes are presented in Supplementary Figure S2.

We used the FRKM (RPKM) method [21] to determine DEGs between compatible and incompatible interactions. We identified 10550 DEGs with a FDR of less than 0.001 and an expression difference of at least two-fold (Figure 3).

We mapped the DEGs to our custom database constructed by merging the pathogen and host genomes. Local Blastx analysis (e-value < 0.00001) revealed that 6072 of the DEGs aligned best to the Pucciniales genome, and 3028 of the DEGs aligned best to the wheat genome. The remaining DEGs yielded no significant homology to either the pathogen or the host; hence, they are considered to be unknown novel DEGs. The 6220 Pst-mapped DEGs (PstDEGs) were selected as the differentially expressed, pathogen-associated candidate genes to proceed with for secretome analysis.

The identified PstDEGs were subjected to Blast2GO analysis to elucidate their functional classification (Supplementary Table S2). A graphical view of the molecular functions and biological processes of the PstDEGs is displayed in Figure 4. We were able to annotate the majority (75.8%) of PstDEGs with their corresponding functional classification. The distribution of molecular function between proteins with catalytic activity (61.73%) and binding activity (61.07%) attributes were practically equal in the Blast2GO level 2 analysis (GO specificity increases as the ‘level’ of the analysis gets higher). During level 3 analysis of the PstDEGs, we observed that the proteins with catalytic activity belonged to the hydrolase, transferase, and oxidoreductase categories (Figure 4A). The Blast2GO analysis illustrates that the pathogen undergoes radical changes during disease progression. A large number of PstDEGs (6072 genes compared with the total predicted protein number of Pst genes which is over 20000) were sorted into various categories of biological processes (Figure 4B). Therefore, we wanted to determine the number of PstDEGs that sorted as enzymes. In total, we cataloged 1848 PstDEGs as enzymes. Among them, hydrolases were the most common, followed by transferases and oxidoreductases (Figure 4C).

We made use of a previously constructed pipeline defined and validated by several research groups [6,23,24]. We extracted the secretome repertoire of 227 SSPs of Puccinia striiformis, which are referred to as PstDESSPs, by selecting all of the DEGs identified in our transcriptome data (6072 DEGs). We identified other SSPs by applying the same pipeline on all of the available genome sequences of Puccinia species: Pst, Pgt, and Ptt (http://www.broadinstitute.org/). The mRNA levels are not directly proportional to the level of protein expressions; nevertheless, transcriptome and proteome expression values are correlated [47]. Moreover, there are reports stating that the overall predicted CSEPs compose a significant proportion of the genome [10,48]. In a similar fashion, we hypothesize that the number of CSEPs during pathogenesis should denote a large amount as well. We then compared all of the identified secretomes in a method similar to that used in the study published by Xia and colleagues, 2017 [10]. Among 227 small secreted proteins, 94 were predicted as effectors using ‘Effector 2.0’ software, and were then classified as Pst-Small Secreted Candidate Effectors (PstSSCEs) (Figure 5A). The percentage of PstDESSPs within the dataset of all DEGs was 3.7% (227 out of 6072). This ratio appeared significantly smaller when compared with the hypothetical prediction of small secreted proteins within PST-78 (6.5%, 1332 out of 20482) and the other Pucciniale genomes. Therefore, we concluded that the number of actively involved small secreted proteins at 10 dpi is indeed much lower than what was hypothetically expected. This outcome indicates that the hypothetical number of PstDESSPs can be exaggerated when obtained using only transcriptome predictions from genome sequences.

Certain types of gene functions are useful for the fungus in pathogen–host interactions, including CAZymes, proteases, peroxidases, and lipases. The PstDESSPs were compared against several reliable public databases (MEROPs, dbCAN, LED, fPoxDB) in order to predict their functions. We identified four proteases and protease inhibitors, four CAZymes, nine lipases, and eight peroxidases among the PstDESSPs (Figure 5B). We searched the Pathogen–Host Interaction (PHI) database to predict the virulence attribute of each PstDESSP based on their sequence similarity. A total of 24 of the PstDESSPs exhibited significant similarity to at least one entry in the PHI database. For example, Unigene31938 (Pstg04010) had a predicted ‘copper-zinc superoxide dismutase’ domain match with PHI383 (C. albicans) and PHI6412 (P. striiformis) (Supplementary Table S3). Although the level of similarity barely exceeded the cut-off value, both these entries have a superoxide dismutase domain. PHI6412 (PsSOD1) has been reported to increase the resistance of the pathogen against host-derived oxidative stress [49]. The results suggest that all of the PstDESSPs that matched with entries in the PHI database are worthy enough of further exploration employing in vivo systems. Hence, a comparison against multiple databases provides a better way to discover CSEPs and understand the importance of the PstDESSPs in Supplementary Table S3. Furthermore, a considerable percentage (80%) of the PstDESSPs remain uncharacterized (Figure 5B), due to the tendency of the candidate effectors to only show homology to the pathogenicity-related domains and not to other typical functional domains.

To understand the relationship between PstDESSPs, we generated a phylogenetic tree using an analysis of multiple sequence alignments. The outcome suggests that there are no significant conserved sequences among the 227 PstDESSPs. In general, we identified three main branches (or classes) of PstDESSPs based on sequence similarities. The purple branch in Figure 5C possesses a substantial amount of the PstDESSPs and is also subdivided into two smaller arms. No significant conserved motifs were detected among the PstDESSPs, except for the Y/F/WxC motif (Supplementary Table S3), which is common to the powdery mildew [50]. The search for the conserved sequences, correlations, and motifs are at the center of the studies using genome-wide analysis or transcriptome-based analysis of CSEPs [51]. Since the discovery of RXLR or LXLFLAK motifs aid the mining of large effector families in oomycetes [52]. Only [YFW]xC motif has been discovered in rust and mildew fungus; however, there is no report for its contribution to the function or localization [48,50,51]. Similarly, we here checked the sequence resemblance and presented it as a phylogenetic tree. Motif mining through MEME motif finder also did not give any significant repeating motif pattern. Hence, we reported the negative result in parallel with the literature.

The PstDESSPs determined in the present study were compared with other published expression datasets of Pst to investigate whether the SSPs in the present study were unique or similar to haustorial differentially expressed effector candidates. Cantu and colleagues studied the expression data for effectors expressed in infected leaf and haustoria samples and provided a pioneer list of SSP tribes in their analysis [7,8]. We analyzed the similarity between tribes of SSPs via the Blastp program (E-value < 0.00001). A total of 116 PstDESSPs were found to share significant homology to haustoria-expressed secreted effector proteins (HESPs), and 132 of PstDESSPs were found to be homologous to infected leaves-expressed secreted effector proteins (IESPs) (Supplementary Table S3). The remaining 94 PstDESSPs that did not match with HESPs or IESPs were found to be differentially expressed at 10 dpi in the present study. Comparison with the list of HSPs published by Garnica and colleagues indicated that 94 of the PstDESSPs are similar to HSPs [9]. Overall, the number of PstDESSPs that we identified that are not similar to those in previously reported datasets is 31, signifying their exclusivity to the present study. Xia and colleagues scrutinized the effector candidates for avirulence (Avr) by comparing the expression datasets of Cantu and colleagues, 2013 [8], [10]. We compared our list of PstDESSPs with the list of Avr candidates, and the results suggest that three previously overlooked Avr candidates and six known Avr candidates are present among our PstDESSP dataset (Supplementary Table S5).

The cysteine content of a candidate effector protein has long been a determinant for making effector predictions [8–10,23,24]. Although several effectors have been discovered to be not rich in their cysteine composition, such as AvrM [53] and AvrL567 [54], we examined the PstDESSPs in terms of cysteine possession. We identified 64 cysteine-rich small-secreted effector proteins (SCEPs) that have four or more cysteine amino acid residues in their mature protein length (not including the N-terminal signal peptide) (Supplementary Table S5).

To find the most promising subset of candidate effectors, we focused on the PstDESSPs that produced positive results in the EffectorP analysis, had cysteine residue counts of more than three, and had homologous matches within the previously mentioned reports (HESPs, IESPs, and HSPs). A total of 14 PstDESSPs met all five criteria and were defined as reliable CSEPs for further investigations (Figure 5D). A more detailed description of the features of the 14 CSEPs is provided in Table 4.

One of the CSEPs from Table 4 was chosen at random to test the reliability of the presented PstDESSPs data. Unigene17495 (Pstg_10917) was predicted to have a transit peptide sequence following the N-terminus secretion signal. Thus, if the effector candidate succeeds in translocation within the host cell, it may be able to partially target the chloroplast, the site of control for programmed cell death (PCD) against both biotic and abiotic stressors. The heterologous expression of Unigene17495 with a GFP tag on N. benthamiana leaves showed chloroplast localization (Figure 6). Its localization to the chloroplast despite being predicted for secretion outside the cell was unexpected. However, re-checking the predictions by other algorithms; ‘TargetP’ in addition to predicting a secretion signal, only after its removal, it also detected a transit peptide within the remainder sequence (111 a.a.). On the other hand, ‘Localizer’ analysis predicted a transit peptide overlapping the secretion signal region. Therefore, the experimentally shown chloroplast localization, instead of originally predicted secretion outside of the cell, can be explained by the presence of a transit peptide and a possible disruption of the secretion signal due to a particular condition in the vicinity of the molecule. Alternatively, in the natural host, trafficking in both locations may be possible or there is a slight change in the physiological condition that may direct a preference of one target location to another.

Unigene17495 (Pstg10917) is conserved among Pst relatives. However, none of the homologs have been previously studied for their function. Interestingly, an effector candidate of the soybean rust pathogen (Phakopsora pachyrhizi), PpEC82 [55], resembles Unigene17495 (Pstg_10917) with an E-value of 1×e⁻¹⁴ (10⁻¹⁴). PpEC82 was reported to localize to the cytosol and nucleus with strong aggregation [56] properties, unlike the chloroplast-targeting Unigene17495 (Pstg_10917) as we have demonstrated. Moreover, Qi and colleagues stated that PpEC82 was able to suppress BAX-induced cell death when expressed in yeast cells [56]. Correspondingly, we conducted cell death suppression assays to test whether Unigene17495 (Pstg_10917) has suppression abilities similar to PpEC82. We used a cell death elicitor (INF1) as a cell death inducer and categorized the resulting suppression of cell death observations into four levels: strong suppression, moderate suppression, weak suppression, or no suppression at all. Across 15 replicates of Unigene17495-treated N. benthamiana leaves, we observed three cases of strong, three cases of moderate, and three cases of weak cell death suppression (Figure 7, Supplementary Figure S3). The other six leaves presented no significant level of cell death suppression. In all 15 trials, leaves treated with the negative controls (GFP and SP-GFP) resulted in tissue collapse and cell death in the infiltrated area. Therefore, we concluded that Unigene17495 (Pstg_10917) without SP, was able to partially suppress programmed cell death triggered by the INF1 elicitor (Figure 7). We need to emphasize that not all attempts produced equal levels of cell death suppression. Although the expression of the effector was verified by microscope analysis at 2 dpi for all of the leaves, 6 out of the 15 leaves revealed no detectable suppression of INF1-induced cell death. Fisher’s exact test calculated a P-value of 0.0007 for the findings, and thus the results were significant (P<0.05). The Unigene17495 (Pstg10917) is up-regulated on susceptible cultivar during the compatible interaction. Our findings suggest that Pstg10917 may be involved in overcoming the host’s defense mechanisms by modulating or suppressing hypersensitive cell death triggered by pathogen recognition. Together, these data show that Pstg10917ΔSP-GFP (Unigene17495) is capable of significantly suppressing PCD induced by INF1.

It should be noted that the complex nature of heterologous expression systems serves as a double-edged sword due to the oscillations in cell death response with N. benthamiana. However, the versatility of the system provides a practical convenience to study candidate effectors. Other candidate effectors of Pst have been studied using heterologous expression systems. For example, the Pst effector PstHa5a23 was reported to suppress PCD triggered by BAX, INF1, and two mitogen-activated protein kinases (MKK1 and NPK1) involved in host cell resistance [57]. PstHa5a23 is a homolog to Pstg_00676 identified in our PstDESSPs list. In a recent publication, Pstg_14695 and Pstg_09266 were able to suppress PCD as part of effector-triggered immunity [58]. Furthermore, Petre and colleagues utilized mass spectroscopy analysis to decipher host cell factors that interact with GFP-tagged candidate effectors in co-immunoprecipitation samples [59]. Eleven of the twenty candidate effectors reported in that study are homologous to our PstDESSPs. These recent reports support the reliability of our PstDESSPs data as a trustworthy collection of candidate effectors that can be used for effector function studies in the future.

The obligate biotrophic fungus Pst causes yellow rust disease by establishing an interface between the host and pathogen using small secreted proteins. Transcriptome sequencing has proven a useful strategy in elucidating the cross-talk of proteins between the host and pathogen, including the interaction between wheat and Pst. We compared compatible and incompatible interaction points in order to identify the DEGs (PstDEGs) at 10 dpi. In silico predictions highlighted the SSPs within the group of PstDEGs. The annotations, comparisons, and characterizations contributed to the identification of a list of the most promising subset of SSPs, that function as phytopathogen effectors. The presented SSP data represent candidates for the functional analysis to understand yellow rust disease in wheat better. In particular, we were able to determine that Unigene17495 (Pstg10917) of our PstDESSPs list is localized to the chloroplasts of infected plant cells and is capable of halting programmed cell death. The remaining 13 promising candidates are to be investigated in future studies to test whether they are indeed effectors of Pst.

All raw data and sequence reads are uploaded on the NCBI BioProject database under the accession number PRJNA545454. The annotations and analysis results are available in Supplementary files.

The authors declare that there are no competing interests associated with the manuscript.

This work was supported by the grant number call FA1208 with the grant number COST-113Z350 to the PI (Mahinur S. Akkaya); and the TUBITAK-BIDEB: [grant numbers 2215, 2211/C (to Bayantes Dagvadorj, Ahmet Caglar Ozketen and Ayse Andac-Ozketen, respectively)].

Ahmet Caglar Ozketen: Methodology, Investigation, Data Curation, Validation, Writing – original draft preparation, Writing – review and editing. Ayse Andac-Ozketen: Investigation, Validation. Bayantes Dagvadorj: Investigation. Burak Demiralay: Data Curation. Mahinur S. Akkaya: Methodology, Funding acquisition, Supervision, Writing – review and editing.

The authors want to thank the Kamoun Lab, Sainsbury Lab, and Tolga Bozkurt of Imperial College London, Department of Life Sciences, London SW7 2AZ, U.K. for vectors used in the study.

BGI

Beijing Genomics Institute

CAZyme

carbohydrate-active enzyme

COG

Clusters of Orthologous Groups

CSEP

candidate-secreted effector protein

DEG

differentially expressed gene/unigene

DESSP

Differentially Expressed Small-Secreted Protein

dpi

days post-inoculation

EST

expressed sequence tag

FDR

false discovery rate

FPKM

fragments per kb per million reads

GO

gene ontology

HESP

haustoria-expressed secreted effector protein

HSP

haustoria-secreted protein

IESP

infection-expressed secreted proteins

KEGG

Kyoto Encyclopedia of Genes and Genomes

NGS

next-generation sequencing

PAMP

pathogen-associated molecular pattern

Pgt

Puccinia graminis f. sp. tritici

PHI

Pathogen–Host Interaction database

Pst

Puccinia striiformis f. sp. tritici

PTI

PAMP-triggered immunity

Ptt

Puccinia triticina

RIN

RNA integrity number

RPKM

reads per kb per million reads

SCEP

small cysteine-rich effector protein

SSP

small-secreted protein