Transcriptome Analysis Based on Lr19-Virulent Mutants Strain Provides Clues for the Pathogenicity-Related Genes and Effectors of Puccinia Triticina

Background: Wheat leaf rust caused by Puccinia triticina (Pt) remains one of the most destructive diseases of common wheat worldwide. Cultivating resistant cultivars is an effective way to control this disease, but race-specic resistance can be overcome quickly due to the rapid evolution of Pt populations. The critical to control wheat leaf rust is to understand the pathogenicity mechanisms of Pt. Results: In this study, the spores of the Pt race PHNT (wheat leaf rust resistance gene Lr19-avirulent isolate) were mutagenized with ethyl methanesulfonate (EMS) and two Pt Lr19-virulent mutants named M1 and M2 were isolated, suggesting that they carry mutations affecting the Lr19-specic avirulent factor. RNA sequencing was performed on samples collected from the wheat cultivars Chinese Spring and TcLr19 that were infected by wild-type (WT) PHNT and the two Lr19-virulent mutant isolates at 14 days post-inoculation (dpi). The assembled transcriptome data were compared to the reference genome “Pt 1-1 BBBD Race 1.” A total of 216 differentially expressed genes (DEGs) were found from three different sample comparisons including M1-vs-WT, M2-vs-WT, and M1-vs-M2. Of these DEGs, 29 common DEGs were shared between M1-vs-WT and M2-vs-WT comparisons. Among the 216 DEGs encoding proteins, 30 were predicted to be effector candidates. Then 6 effector candidates (PTTG_27844, PTTG_05290, PTTG_27401, PTTG_27224, PTTG_26282, PTTG_25521) were veried that these genes were differentially expressed during Pt infection by quantitative real-time PCR (qRT-PCR) and were validated on tobacco, and the results showed that PTTG_27401 could inhibit progress of cell death (PCD) induced by BAX. Conclusions: Our results showed that a large number of genes participate in the interaction between Pt and TcLr19, which will provide valuable resources for the identication and targeting of AvrLr19 effector candidates and pathogenesis-related genes.


Background
Wheat leaf rust caused by Puccinia triticina (Pt) is one of the most common and severe diseases in the wheatgrowing regions worldwide. The yield losses of wheat infected with Pt ranges from 5% − 15% and the yield of wheat infected with leaf rust at the seedling stage can be reduced by 50% or even more [1]. Considering the impacts of global climate change, temperature and humidity conditions may be more suitable for the proliferation and epidemic of Pt in the future which will further deleteriously impact wheat yields. Control of leaf rust mainly relies on the fungicide application and deployment of cultivars carrying resistance genes. Genetic resistance is the most effective, environmentally safe, and economically feasible approach to reduce the damage caused by Pt. However, monoculture of select resistant varieties leads to the host selection pressures that drive Pt evolution and promote the continuous emergence of new toxic races, which often leads to the decline of wheat resistant varieties after several years of planting.
To date, over 80 genes conferring leaf rust resistance derived from Triticum aestivum wheat cultivars, wild wheat, grass species, and durum wheat have been identi ed and designated as Lr1 to Lr79 [2]. Among them, only Lr9, Lr19, Lr24, Lr34, Lr37, Lr38, Lr46, Lr47, Lr51, Lr53, and Lr68 confer effective resistance to leaf rust currently. Lr19 is derived from the grass Agropyron elongatum and was transferred to the long arm of wheat chromosome 7D [3]. The wheat cultivar carrying Lr19 is still effective against all Pt races in Asia, Australia, Canada, and Europe, and resistance is expressed during the whole growth period [4,5]. Therefore, Lr19 holds the potential to be deployed in combination with other Lr genes in the eld to confer durable resistance against leaf rust worldwide [4,6]. It has been reported that Lr19 is associated with increased grain yield, which promotes it as an important gene for wheat breeding against Pt-mediate yield losses [7]. To prevent Lr19 from being out-competed by newly emergent virulent races, cloning avirulent gene AvrLr19 that can be recognized by wheat leaf rust resistance gene Lr19, determining the molecular mechanism of AvrLr19 and Lr19 interactions, and monitoring the natural Pt population changes in response to AvrLr19-Lr19 resistance are key to enable the long-term deployment of Lr19 in leaf rustresistant wheat varieties. At present, there is no virulent strain of Lr19 in the eld, which not only bene ts the use of Lr19 but also restricts the development of AvrLr19.
Mutation is the most important avenue for creating new rust races and genotypes. Ethyl methanesulfonate (EMS) is an alkylating chemical mutagen that generates single nucleotide polymorphisms (SNPs) and insertions and deletions (Indels), resulting in amino acid sequence variation, and nally leads to phenotype changes that can be selected upon. Mutagenesis integrated with genomic sequencing is an e cient way to study the relationships between phenotypic traits and associated genotypes, leading to the identi cation of fungal effectors or Avr genes. For example, Salcedo et al. obtained stripe rust mutants through EMS-induced mutation and performed genome sequencing to obtain AvrSr35 candidate genes, and then veri ed the candidate genes through co-expression of AvrSr35 and its corresponding resistance gene in tobacco and wheat to trigger cell death [8]. Li et al. screened 30 mutant variants from the least virulent isolate generated by EMS mutagenesis and candidates Avr genes were determined by sequencing [9]. Transcriptomics has proven to be an instrumental molecular tool to help identify virulence effectors and Avr genes [10][11][12]. As the host and pathogen interact in a battle for supremacy, the underlying transcriptional regulation of gene expression in the plant and pathogen provides clues to their defense and virulence mechanisms, respectively [13].
The gene to gene hypothesis proposed by Flor indicates that only the host with an "R gene" is resistant to the homologous Avr gene in pathogens [14]. Jone proposed the famous "ZigZag" model in 2006 to analyze the molecular mechanism of interaction between plants and their pathogens [15]. In order to inhibit plant defense responses, pathogens secrete a series of effector proteins through the haustorium to interfere with or ablate the plant defense response, so as to meet their own growth needs. In recent years, with the continuous improvement of sequencing technology and the reduction of sequencing costs, and the development and application of prediction software such as SigalP [16], TargetP [17], TMHMM [18], PredGPI [19] and Pfam [20] have improved the screening e ciency of candidate Avr genes of rust. To date, a handful of Avr genes have been identi ed in rust pathogens, including AvrL567, AvrP123, AvrP4, AvrM, AvrL2, and AvrM14 from the ax rust pathogen [21], as well as RTP1 in bean rust [22], and PGTAUSPE10-1 from wheat stem rust [23]. At the end of 2017, two articles reported the successful cloning of the Avr genes AvrSr35 [8] and AvrSr50 [24] of wheat stem rust, and preliminarily revealed their interaction with corresponding resistance genes. Some effector candidates for Lr26, Lr9, and Lr24 materials were obtained from Pt [25] and 20 effector candidates of Lr20 [26], but their biological functions have not been determined, so, no known Avr genes have been identi ed in Pt so far.
In this study, we aim to use an EMS mutagenized Pt race PHNT to obtain Lr19-virulent mutants, and identify differentially expressed genes (DEGs) associated with the Pt infection. Our results provide resources for identi cation of AvrLr19 candidates and pathogenicity-related genes, and lay a foundation for revealing the pathogenic mechanism of Pt.

Results
Lr19 triggers a resistance response at early stages of infection Leaf rust resistance conferred by Lr19 is best expressed in all stages of the wheat plant and culminates in a hypersensitive response (HR) at the infection site which is also known as race-speci c resistance [5]. In order to con rm the recognition between Lr19 and its respective Avr gene in the Lr19-avirulent Pt race PHNT, the phenotype and histopathology examination of Pt-infected leaf tissues from resistant TcLr19 (Lr19+) and susceptible Chinese Spring (Lr19-) wheat lines were analyzed. Compared to Chinese Spring (Lr19-) wheat lines, the development of fungal infection hyphae (stained blue) stopped before the formation of a haustorium, the structure with which the fungus extracts nutrients from its host plant, in the resistant wheat line TcLr19 (Lr19+).
Imaging at 2 dpi revealed fungal growth in susceptible Chines Spring but no further fungal growth in TcLr19, the dead host cells in TcLr19 were stained green, and no dead cells were revealed in Chinese Spring (Fig. 1). This early immune response is consistent with pronounced HR symptoms and suggested early expression of a fungal gene recognized by Lr19, which demonstrated that Lr19 triggers a resistance response at the early stages of infection.
Lr19-virulent mutants were obtained and con rmed The Lr19-avirulent Pt race PHNT was used to create Lr19-virulent mutants via EMS mutagenesis. The germination rate of PHNT spores under different EMS concentrations was calculated, and the results showed that the germination rate of spores was 50% at an EMS concentration of 0.005 M ( Fig. 2A). Two Pt mutants virulent to the Lr19 allele were isolated, named M1 and M2, suggesting that they carried Lr19-speci c Avr factor (Fig. 2B). Confocal microscopy of wheat leaves from the compatible Chinese Spring cultivar infected with the wild-type (WT) and mutant Pt strains was performed to investigate the effect of the AvrLr19 gene mutations on pathogen virulence. In order to ensure the reliability of experimental materials of TcLr19, molecular markers of Lr19 [27] were used to detect the existence of Lr19 (Fig. 2C). To validate that only Lr19-virulence was affected by the EMS mutagenesis, M1 and M2 were analyzed by infecting a standard full set of near-isogenic lines that carry different Lr genes, with the WT PHNT race inoculated as controls (Fig. 2D, Supplementary Table S1). This panel is broadly used for quick pathotyping of unknown Pt isolates using a standard 4-letter code. Only TcLr19 inoculated with M1 and M2 demonstrated a change from resistant ('0') to susceptible ('3 +'), while the susceptibility of other nearisogenic lines remained unchanged, indicating that AvrLr19 was altered in M1 and M2 alleles and that Lr19virulent mutant strains were successfully obtained.

RNA-seq analyses display the information of Pt race PHNT and mutants
To characterize the DEGs between Chinese Spring-Pt and TcLr19-Pt mutants during infection, RNA-seq of TcLr19 inoculated with different Pt strains were performed. Following inoculation and incubation, leaves with obvious phenotypes were harvested, then the total RNA was extracted from M1, M2, and WT (Fig. 3A) and was used for RNA-seq cDNA library construction. Quality statistics of sequencing data and the percentages of reads for each sample were mapped to the Pt BBBD race reference genome (Supplementary Table S2 and Table S3). The percentage of reads that aligned to the Pt reference in M1, M2, and WT are shown in Fig. 3B, respectively. The proportion of reads aligning to the wheat reference Chinese Spring was higher in WT than that in the mutants  (Table 1). Table 1 Results of KEGG analysis of DEGs between M1-vs-WT and M2-vs-WT.

Class
No   The expression pro les of selected effector candidates were validated by qRT-PCR assay We searched the data of Pt race PHTT(P) [28] and foundthat 6 of 30 effector candidates were differently expressed (Fig. 7A). So they were veri ed the expression levels of at different time post during Pt race PHNT infection by qRT-PCR analysis. The expression level of PTTG_27844 at 0.5 days post-inoculation (dpi) was 3-fold higher than that of germ tubes (GT), 10-fold higher at 1 dpi, and increased to 304-fold higher at 4 dpi. The expression level of PTTG_05290 at 0.5 dpi was 36-fold higher than that of GT, 28-fold higher at 1 dpi, and increased to 218-fold higher at 4 dpi. In addition, the expression patterns of PTTG_27401 and PTTG_25521 were similar to PTTG_27844 and PTTG_05290, but the expression levels of PTTG_27844 and PTTG_25521 were lower. PTTG_26282 expression peaked at 0.5 dpi and 4 dpi. PTTG_27224 expression peaked at 0.5 dpi (Fig. 7B). These results were similar with the gene expression patterns predicted based on their FPKM values [28].  Figure S4). However, there was no necrosis at the sites where BAX and PTTG_27401 were in ltrated together, indicating that PTTG_27401 inhibited BAX action, preventing PCD from being induced, with potential toxic function (Fig. 8).

Discussion
It is well-known that mutation is the ultimate source of genetic variation, providing new alleles and genotypes [29].
The innovation of the present study is that we developed two Lr19-virulent mutants via EMS mutagenesis, M1 and M2, and con rmed that the mutation site only altered the virulence of Pt race PHNT to Lr19, but not other Lr genes, which bene ts the screening of AvrLr19 effector candidates.
RNA-seq has been extensively used to characterize transcriptional changes in both the host and pathogen at different stages of pathogen infection. Thus, it is possible to utilize transcriptomic data to identify DEGs during compatible and incompatible interactions. In this study, 216 DEGs including 29 DEGs found in both M1 and M2 were obtained, which suggested that they may play a certain role in infection of PHNT on TcLr19. PTTG_03515 is annotated as tyrosinase which can form melanin through a series of catalytic reactions. Melanin has been reported to play a role in the pathogenicity of fungi, and it is believed to contribute to fungal structure formation during infection, or to in uence the host response to infection [30]. Another putative protein, PTTG_05196 is annotated as zinc nger protein that plays a role in pathogenic bacteria infection, and it has been proven that it can promote the occurrence of diseases [31].
30 effector candidates with certain criteria were identi ed. These genes provide valuable resources for screening AvrLr19 candidates and pathogenesis-related genes. Among them, PTTG-01075 and PTTG-11895 are annotated as clan D glycoside hydrolases (GH-D), which are commonly referred to as plant cell wall degrading enzymes (PCWDEs) [32]. PTTG-12205 is annotated as a Major Facilitator Superfamily (MFS) of sugar transporters, which facilitates the transport of diverse molecules like sugars, vitamins, amino acids, hormones, etc. across cell membranes. MFS proteins might be relevant to the uptake of nutrients following the lysis of the host plant tissues, as well as the transport of toxins and antimicrobial compounds [33]. PTTG_25982 was not only a common DEG in M1-vs-WT and M2-vs-WT comparisons, but was also predicted to be an effector candidate. Zhao et al sequenced Pt race PHTT(P) and analyzed the functional factors with conserved motifs [28], in this study, four common genes (PTTG_25271, PTTG_04128, PTTG_08504, and PTTG_08503) encoding PNPi like effectors were found.
Salcedo et al. [8] successfully cloned AvrSr35 and analyzed the expression patterns in wheat infected by wheat stem rust at 0, 1, 2 and 4 dpi, the results showed that AvrSr35 expression peaked at 4 dpi. In this study, we identi ed that PTTG_27401 suppressed PCD triggered by BAX and its expression peaked at 4 dpi. Further investigations on the molecular receptor of effector candidates should greatly improve our understanding of the avirulent mechanism of the Pt.
We pro led genes associated with the infection process of the Lr19-virulent mutants Pt strains derived from PHNT. By analyzing the transcriptomes of all samples, we found that the ratio of M1_sample 1, M1-sample 3, and M2-sample 6 to the reference genome were lower than other samples. It may be due to insu cient spores in the sample. It was reported that 15 Sr35-virulent mutants generated by EMS were used to clone AvrSr35 [8]. In the future, more Lr19-virulent mutants may need to be obtained in order to improve the power of this study. In addition, the function of AvrLr19 candidates and pathogenesis-related genes of Pt will be identi ed by using the bacteria type III secretion system and host-induced gene silencing. By de ning AvrLr19 and its key functions, AvrLr19 can be used to develop gene-speci c molecular markers for directly monitoring Pt population changes in the eld, which can be used for real-time monitoring of leaf rust races combined with marker-assisted selection (MAS) of wheat leaf rust resistance genes to predict the resistance of varieties and inform the correct use of resistant varieties. Our results provide valuable resources for identifying and characterizing AvrLr19 effector candidates and pathogenic genes, and lay a foundation for further elucidating the pathogenic mechanism of Pt and analyzing the disease-resistance mechanism of Lr19.

Conclusion
In the present study, we obtained 2 Lr19-virulant mutant isolates developed by EMS mutagenesis of Pt race PHNT. By RNA-seq the WT isolate and the mutants M1 and M2, we aligned the sequences of the mutants with the wildtype isolate. We pro led genes associated with the infection process of the Lr19-virulent mutants Pt strains derived from PHNT. A total of 216 DEGs were found from M1-vs-WT, M2-vs-WT, and M1-vs-M2. Of these DEGs, 29 common DEGs were shared between M1-vs-WT and M2-vs-WT comparisons. Among the 216 DEGs encoding proteins, 30 were predicted to be effector candidates. Among them, 6 effector candidates were chosen and validated expression pattern and validated on tobacco, the results showed that PTTG_27401 could inhibit necrosis induced by BAX. Our results provide valuable resources for identifying and characterizing AvrLr19 effector candidates and pathogenic genes, and lay a foundation for further elucidating the pathogenic mechanism of Pt and analyzing the disease-resistance mechanism of Lr19.

Plant materials and Pt isolates
Wheat seedlings of cultivar TcLr19 (Tc*6/RL6040), which is a near-isogenic line containing the Lr19 gene, were requested from the Cereal Disease Lab of the USDA located at the University of Minnesota and were preserved in the Pt laboratory at Hebei Agricultural University. Wheat seedlings of cultivar Chinese Spring collected from our laboratory. Seedling plants of the all wheat cultivar were grown in a glasshouse. Pt pathotype 07-10-426-1 (PHNT race), collected from China and preserved in our laboratory, was used to inoculate wheat according to Roelfs' standards [34]. Seedlings (14 days old) from wheat lines Chinese Spring (Lr19-) inoculated by spraying with a suspension of Pt urediniospores were used as controls.

Staining and microscopic observation
To identify fungal tissues and the presence of dead cells in the leaf mesophyll of infected plants, the samples were stained via the Rohringer uorescent staining method [35]. Stained tissues were observed under optical microscopy, and at least 10 fungal infection sites from three biological replicates were analyzed for each stained leaf sample.

EMS mutagenesis of Pt spores of the PHNT race
The Lr19-avirulent Pt race PHNT was used to create Lr19-virulent mutants by treatment with EMS as previously described [8]. Four different concentrations of EMS were used for mutagenesis: 0.1 M, 0.05 M, 0.01 M, and 0.005 M. After treated 2 hours with EMS, the germination rate of spores was observed using microscopy and unmodi ed Pt spores were used as control. Spores with at least 50% germination rate were inoculated on 12-dayold wheat TcLr19 seedlings. Inoculated seedlings were incubated in a dew chamber (22 °C) in dark conditions for 12 hours and then were kept in a growth chamber at 22 °C with a 16 h/8 h (day/ night) photoperiod.

Con rmation of mutation site
To ensure that the obtained Lr19-virulent mutant strains were not escapes of other rust strains and to validate that only Lr19-virulence was affected by EMS mutagenesis, the mutants and wild type leaf rust were analyzed by infecting a standard set of differential wheat near-isogenic lines that carry different Lr genes. Wheat seedlings of near-isogenic lines were requested from the Cereal Disease Lab of the USDA located at the University of Minnesota and were preserved in the Pt laboratory at Hebei Agricultural University. After 14 days, the phenotypes of hosts infected with wild-type and mutant strains were recorded according to Roelfs' 6-level (0, ;, 1, 2, 3, 4) classi cation method [34]. In addition, the molecular marker of Lr19 was used to detect the validity of TcLr19 [27].

Transcriptome sequencing
Before sequencing, the Pt culture was puri ed by single-pustule isolation and pathotyped according to Roelfs' standards [34]. 14-day-old seedlings of TcLr19 inoculated with the urediniospores of the PHNT mutants were collected, and susceptible Chinese Spring wheat cultivar inoculated with the urediniospores of Pt race PHNT were used as a control. Plant samples were collected from severely diseased leaves at 14 dpi and three biological replicates were prepared for each group. Transcriptome sequencing and analyses were conducted by OE Biotech Co., Ltd. (Shanghai, China). Total RNA was extracted from the different samples and one microgram of highquality total RNA was taken for RNA-seq library construction. The quantity and quality of extracted RNA were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, UK). Sequencing was carried out on the Illumina HiSeqTM 2500 platform. Raw reads were processed using Trimmomatic [36] to obtain the clean reads.
Then, the clean reads were mapped to reference genome "Pt 1-1 BBBD Race 1" [37] using Hisat2 [38] to obtain the location information on the reference genome or gene, as well as the speci c sequence feature information of the sequenced sample. The FPKM values for each of the extracted transcripts were calculated using Cu inks [39]. The read counts of each gene were obtained by HTSeq-count [40]. DEGs were identi ed using the DESeq (2012) R package [41]. "FDR-adjusted p-value < 0.05" and "|fold Change| >2" were set as the thresholds for signi cantly differential expression. Hierarchical clustering analysis of DEGs was performed to explore gene expression patterns. GO enrichment and KEGG [42] pathway enrichment analysis of DEGs were respectively performed using R based on the hypergeometric distribution.

Effector candidates among the 216 DEGs
To search for Pt effector candidates, we adopted the pipeline that shares common features with the effector prediction pipelines described for lamentous plant pathogens [43]. We rst detected the presence of an Nterminal signal peptide through SignalP 5.0 [16], which indicates that candidates are presumed to have extracellular functions. Next, proteins with predicted sequences containing transmembrane helices were excluded from the set of predicted extracellular proteins using TMHMM 2.0 [18]. Subcellular localization signals was anlyzed by TargetP 2.0 [17]. The presence of a glycosylphosphatidylinositol (GPI) anchor was predicted using PredGPI [19] to identify proteins anchored to the membrane. Conserved protein domains were identi ed using Pfam software [20]. Perl software was used to search for known motifs, including YxSL[R/K] [44] and RxLR [45] motif commonly detected in oomycetes,     The numbers of DEGs are noted in each section of the Venn diagrams.

Figure 5
Bar graph showing the number of 216 DEGs identi ed utilizing gene ontology (GO) enrichment analysis for involvement in speci c molecular function, cellular component and biological processes. In the gure, red color represents molecular function, green color represents molecular function, blue color represents molecular function, the Y coordinate is the name of GO term, and the X coordinate is number of gene.  Transcriptional pro le of genes during the Pt pathotypes infection measured by qRT-PCR. The transcript levels of the six selected effector candidates during Pt infection at 0.5, 1, 4, 6, 8, 11, and 15 dpi were determined by qRT-PCR assay. Samples collected from the GT of Pt uredospores served as a control. The transcript levels for all genes were expressed as linearized fold-PtActin levels, which were calculated according to the formula 2△CT. Data were expressed as mean values ± SE from 3 biological replicates. An asterisk (*P < 0.05,**P < 0.01, ***P <0.001) indicates a signi cant difference between the control and infection samples by Dunnett' s test.