Puccinia triticina pathotypes-THTT and THTS display complex transcript profiles on Thatcher

Background: Wheat leaf rust is an important disease worldwide. Understanding the Puccinia triticina f.sp tritici Pt ) pathogenic molecular mechanism and inconstant toxic region are highly important for managing the disease. The present study aimed to analyze pathogenic divergences between Pt isolates. Results: Total RNA was extracted from Thatcher infected by two Pt isolates, Tc361_1 (THTT) and Tc284_2 (THTS) at 144 hours post inoculation (hpi). The mRNA then was harvested and sequenced. Results indicated that a total of 2,784 differential expressed genes (DEGs) were detected. Forty-five genes were specifically expressed in THTT including transcription initiation factor, transmembrane transporter activity and so on, while 26 genes were specially expressed in THTS including GTPase activity, ABC transporter and so on. Fifty-four differential expressed candidate effectors were screened from the tow isolates. Two candidate effectors were chosen and validated on tabacco, and the result showed that they can inhibit the necrosis induced by Bax. qRT-PCR of 12 significantly different expressed genes were carried out to prove that the results of RNA-seq is similar to that of qRT-PCR at 144 hpi and show their expression level in the early stage and the difference between two Pt pathotypes. Conclusion: The results obtained in this study showed that the transcript profiles between THTT or THTS and Thatcher were very complicated. metabolism, tight junction, viral myocarditis, lysosome, thyroid cancer, fatty acid biosynthesis, other glycan degradation, pathways in cancer, Amino sugar and nucleotide sugar metabolism, MAPK signaling pathway - yeast.

. The haustorium affects the metabolism of the host cell [4,5]. In addition, the haustoria secrete toxic effector molecules into the extra-haustorial matrix which are then transported to the host cell to alter the plant cellular defence, architecture and metabolism, ultimately leading to a compatible plant-pathogen interaction [6,7].
The rust genome was estimated to be 100-135 Mbp [8,9]. Since it is large and 43% of the sequences are repeats, assembly and analysis are challenging. Thara et al. [10] used suppression subtractive hybridization (SSH) to identify 69 genes induced in the fungus as it developed in plants four days after inoculation. The proteome of the susceptible wheat/Pt interaction was interrogated using twodimensional polyacrylamide gel electrophoresis (2D-PAGE) at nine days post-inoculation [11]. Hu et al. [12] developed an expressed sequence tag (EST) database, which represent each of several lifecycle stages of Pt. Panwar et al. [13] found that the PtMAPK was involved in pathogen development, as are its orthologues in Ustilago maydis. Segovia [14] found 10 potential avirulence candidates in 432 EST's which were derived from haustoria and infected plants, when she identified the wheat leaf rust genes expressed during the early stages if infection.
RNA-seq is beneficial for revealing information regarding the transcriptome of Pt. Duplessis et al. [15] sequenced Mlp 98AG31 of M. larici-populina and Pgt SCCL and found 16,399 and 17,773 secreted proteins, respectively, and the annotated functions of these proteins were identified. Cantu et al. [16] identified five candidate effector proteins with polymorphisms in 2,999 secreted proteins by sequencing the different virulence of stripe rust races PST-87/7 and PST-08/21 from the UK. Bruce et al. [17] inoculated wheat with six different Pt races and then sequenced wheat leaves severely infected on six days after inoculation using the Illumina platform. Five hundred thirty-two candidate secreted proteins were predicated, of which 456 proteins were present in all the tests of races and twelve candidate avirulent genes were predicated. Recently, a comparative genomic approach was integrated with an association analysis to identify candidate effector genes corresponding to Lr20 in phenotype-paired Pt isolates from Australia [18]. Twenty Pt isolates from Australia, comprising 10 phenotype-matched pairs contrasting Lr20 pathogenicity were analyzed using whole genome sequencing. This was the first report to integrate phenotype-genotype associations with effector prediction in Pt genomes which is an approach that may circumvent technical difficulties in working with obligate Pt and accelerate the identification of avirulent genes. The AvrSr35 and AvrSr50 were cloned using the mutation, the RNA-seq and genome sequencing of the isolates of Puccinia graminis f.sp tritici (Pgt) [19,20]. The discovery of AvrSr35 and AvrSr50 not only provides a new tool for the surveillance of Pgt identification of host susceptibility targets and characterization of the molecular determinants of immunity in wheat but also indicated the complication of pathogenic mechanism of Pgt. Relevant studies on investigating candidate effector proteins have been performed internationally and have a great progress on rust fungi. But there is limited understanding regarding the pathogenic mechanism. Therefore, understanding the differential expression of different Pt races on the same susceptible plant is important for revealing the difference on pathogenecity and pathogenicity differences of Chinese Pt isolates. To gain more insight into the molecular basis of Pt/wheat biotrophism and to identify genes involved in pathogenicity and virulence, RNA-seq were conducted to analyze the differentially expressed genes between two Pt pathotypes THTT and THTS at 144 hpi after interacting with Thatcher.

RNA preparation and establishment of transcriptome library
A pot with equal length in diameter and height of 10 cm was used to plant approximately 20 seeds of the susceptible wheat cultivar, Zhengzhou 5389 (provided by our laboratory) in the greenhouse.
Seven-days old seedlings were artificially inoculated separately with Pt pathotypes-THTT and THTS (provided by our laboratory), which are epidemic strains with only differences between "T" and "S".
( Table 1). The inoculated plants were placed in a moist chamber overnight at 18-20°C in dark and then transferred to greenhouse chambers and cultured at 20±2°C for 12 hours' light. The Pt pathotypes were purified, propagated in Zhengzhou 5389 and then inoculated in the susceptible Thatcher line. We sampled the same position of the inoculated Thatcher leaf tissue at 144 hpi for total RNA of Pt respectively. The RNA was isolated using an RNA isolation kit (RNeasy Plant Mini Kit produced in Tiangen, Beijing) according to the manufacturer's instructions. The Illumina HiSeq 2000 sequencing platform was used for sequencing total RNA (sample IDs areTc361_1 and Tc284_2 respectively) after testing the quality of the extracted RNA (conducted by BGI Company).

Data analysis
The differentially expressed genes (DEGs) according to the sequencing data were screened using the RPKM algorithm [21] (Reads Per Kb per Million reads) and threshold value. We considered genes with FDRs less than or equal to 0.001 and the difference is less than 2-fold as DEGs, along with their annotated function and signalling pathway according to Gene Ontology (GO) and KEGG Library analysis, were prepared for the cluster analysis. The secreted proteins were predicted by SignalP 4.1,

qRT-PCR analysis
The expression of the differentially expressed genes and their primers for qRT-PCR are presented in Table S1. The leaves were infected with the two different pathogenic types of Pt, THTT and THTS respectively, and sampled at 0 hpi, 6 hpi, 12 hpi, 18

Agrobacterium-mediated transient transformation of Nicotiana benthamiana
Two candidate effectors were chosen to identify if they could inhibit the cell death induced by Bax.
The constructed recombinant plasmid PGR107::Unigene17565 and PGR107::Unigene23118 using primers extended with a ClaI and an XmaI site, respectively (Table S2), was transferred into Agrobacterium competent cell GV3101, and positive clones were picked and shake at 28 ° C. Diluting them and the cell death inducer (Bax) to OD 600 = 0.3. The 2nd to 4th leaves of the 6-8 week-old N.
benthamiana plant leaves were infiltrated by effectors and Bax at the same time. Leaves were photographed and decolored in 7 days after inoculation.

Plant response to infection
The differences found between the strain THTS and strain THTT included only "S" and "T" in the corresponding pathogenic types which means THTS was avirulent to Lr18 in the virulence formula, while THTT was virulent to Lr18. However, THTT was avirulent to Lr36 and Lr44 in the virulence formula, while THTS was virulent to them (Table 1).

Differentially expressed genes based on RNA-seq
The total reads from the data base for THTS and THTT were 7,206,771 and 7,040,346. The gene expression levels were compared between THTT as the treatment group and THTS as the control group. Twenty-one hundred seventy-two genes with different expression levels including 1,367 genes that were specifically expressed in THTT and 1,222 genes that were specifically expressed in THTS were screened. And 2,784 significantly DEGs (FDR ≤0.001 and a difference ≥2 times), of which 1,708 were up-regulated and 1,076 were down-regulated ( Figure 1A).were obtained. And there are 45 and 26 DEGs expressed specifically in THTT and THTS, respectively. ( Figure 1B). Few of them were annotated including zinc metalloprotease, phospholipase, RNA polymerase, cell wall protein, actin cytoskeleton-regulatory complex protein and glycoprotein glucosyltransferase among up-regulated genes (Table S3), and phosphate acyltransferase, dihydroneopterin aldolase and ATPase among down-regulated genes (Table S4).
In THTT VS THTS, 2,784 DEGs were enriched in the three types of Ontology functional classifications, including biological processes, cellular components and molecular functions ( Figure 2A and Table 2).
A total of 875 genes were annotated to biological process, 539 genes were annotated to cellular component and 971 genes were annotated to molecular function.
A Q value less than or equal to 0.05 indicates that the DEGs were significantly enriched in the KEGG pathway. In total, 1,562 genes of 2,784 DEGs were annotated to 153 KEGG pathways. Pathways with higher confidence in the enrichment were shown in Figure 2B including ribosome, starch and sucrose metabolism, tight junction, viral myocarditis, lysosome, thyroid cancer, fatty acid biosynthesis, other glycan degradation, pathways in cancer, Amino sugar and nucleotide sugar metabolism, MAPK signaling pathway -yeast.

Candidate effectors in 2,784 DEGs
One hundred fifty-nine proteins containing signal peptides were predicted among 2,784 proteins coded by DEGs using Signal P 4.0. A total of 118 proteins were found by TMHMM 2.0 without transmembrane domain. Then Effector P 2.0 [22] is used to further clarify and 54 candidate effectors were screened. Compared to THTS, there were 21 up-regulated candidate effectors genes and 33 down-regulated candidate effectors genes. The structural analysis was conducted by bioinformatics method (Table S5). Their size range from 59 to 237aa. Cysteine content of 41 candidate effectors are more than three. The predicted 54 effector proteins were analyzed using the Nr database. The results showed that only Unigene22186 had functional annotation, which was copper/zinc superoxide dismutase, and the rest were all hypothetical proteins. The domain analysis was performed using

qRT-PCR analysis of 12 differentially expressed genes
The result of RNA-seq is similar to that of qRT-PCR at 144 hpi. The expression trends of genes are related to the formation of structures after rust infection. By 12 hpi, germ tube, appressorium and substomatal vesicle were clearly formed. Primary hyphae and haustorial mother cells were formed at 18-24 hpi, secondary hyphae and haustora were formed at 36 hpi. At 48 hpi, small mycelial mass was formed [23]. In order to confirm that the trend of gene expression is consistent with the result of RNAseq and discover the expression characters at different time points, we selected differentially expressed important genes for real-time fluorescence quantitative test ( Table 2).
Three genes enriched in different KEGG pathways were chosen to analyze. CL622.Contig1 ( Figure 3A) in THTT peaked at 288 hpi and, in THTS, peaked at 144 hpi. CL3900.Contig1 ( Figure 3B) was significantly up-regulated at 6 hpi in THTT and THTS, with different expression quantity which was significantly higher in THTT than in THTS, followed by a decrease. This gene was annotated to adenosinetriphosphatase. Based on the expression properties, the spores were in the stage of the formation of germination and bud tube at 6 hpi, so we speculate that the CL3900.Contig1 was thought to act as energy carriers that participate in the germination or the formation of germ tubes.
Unigene18070 ( Figure 3C) was annotated to be involved in helicase activity and DNA binding.
According to the qRT-PCR results, in THTS it reached the expression top level at 24 hpi while in THTT at 12 hpi and the expression quantity was higher than in THTT. Huang et al. [23] observed that Pt had formed appressoria since 8 hpi. Additionally, haustorial mother cells and primary hyphae developed when Pt contacts an epidermal or a mesophyll cell [24]. We suggested that the appressoria formation was earlier in THTT than that in THTS.  Figure 3I with superoxide dismutase activity were similar in both THTS and THTT and the expression level in THTS was also higher than that in THTT.
Two peaks were observed at 12 hpi and 144 hpi in THTS, the highest peak appeared at 144 hpi. In THTT, the expression began to increase at 48 hpi and reached its peak at 144 hpi. The expression levels of CL6956.Contig1 ( Figure 3J) annotated as hydrolase activity, was higher in THTT than those in THTS. These genes all peaked at 12 hpi and then gradually decreased.
The other two differentially expressed effectors were also analyzed. The expression of Unigene11935 ( Figure 3K) is very low in THTT at every time point, but peaked at 36 hpi in THTS. There was no expression of Unigene11683 ( Figure 3L) in THTT at 144 hpi. However, two peaks were detected at 12 hpi and 36 hpi. And reached its peak at 48 hpi in THTS. Therefore, Unigene11935, which only expressed in THTS, is a specially DEG.

Bax onN. benthamiana
To test whether the effector proteins can inhibit PCD, the combined plasmids, Bax, negative control-eGFP were injected into N. benthamiana. It is clearly that there was no HR at the site where the vector, Unigene23118 and Unigene17565 were injected alone, respectively. And the necrosis occurred in the sites where Bax alone and Bax and vector together were injected. However, there was no necrosis at the site of Bax and whether Unigene23118 ( Figure 4A) or Unigene17565 ( Figure 4B), indicating that the gene inhibited the action of Bax and PCD could not be produced. We can come a conclusion that Unigene23118 and Unigene17565 have the function of effectors.

THTT and THTS expressed complicated expression patterns when interacting with Thatcher at 144hpi
Pathotypes expressed in many differences in the interaction process with Thatcher when we detected at 144 hpi, including 1,076 genes that were significantly up-regulated of the 10,483 up-regulated genes and 1,708 genes that were significantly down-regulated of the 10,689 down-regulated genes in (THTT) compared to (THTS). We aligned 566 significantly up-regulated genes and 1,225 down-regulated genes to KEGG and GO. Of these 2,784 genes with significantly different expression, 100 genes were present in the Nr NCBI Library, and 78 genes of them were significantly up-regulated, while 22 genes were significantly down-regulated in THTT compared to THTS. A GO analysis was performed to classify these differentially expressed genes into 3 major biological categories. These genes belonged to 19 classes involved in biological processes, 9 classes involved in cell location categories and 12 classes involved in molecular functions.
According to the RNA-seq library, CL622.Contig3 was annotated as mannosyl-oligosaccharide glucosidase (OsMOGS). OsMOGS was involved in N-glycan synthesis in rice, which played an important role in establishing and maintaining the growth of auxin and the growth and development of the root system [25]. Therefore, we speculated that this gene participated in the development of epithelial cells and mycelium growth.
Unigene1676 is annotated to ubiquitin E3 ligase. E3 ubiquitin ligases transfer Ub to one or more Lys residues in the substrate by linking the C-terminal Gly of Ub with a Lys of the target protein (and/or a Lys of the Ub itself). Ubiquitin E3 ligase is closely related to the regulation of the cell cycle, tumourigenesis, cell proliferation, cell apoptosis, signal transduction, cell growth, cell immunity, inflammation and the regulation of the replication and repair of DNA. F-box protein (SCF) E3 ligases are the largest E3 gene family, of which the F-box protein is the key component to determine substrate specificity. Some studies have revealed that the deletion mutant of GrrA, a F-box protein in Aspergillus nidulans is unable to produce mature ascospores because of a block in meiosis [26].
GTPases (singular GTPase) are a large family of hydrolaseenzymes that can bind and hydrolyze GTP.
The GTP binding and hydrolysis takes place in the highly conserved G domain common to all GTPases.
Zhang et al. [27] verified the functions of all six Rho GTPases in Fusarium graminearum by constructed the deletion vectors. Mutation △Fgrac1, △Fgcdc42 and △Fgrho4 were drastically reduced in growth rate and △Fgrho4 lacked aerial hyphae. △Fgrho2 and △Fgrho3 were less drastic on CM plates compared to other mutants. FgRho1 is essential for fungal survival. FgRho2, FgRho4, FgCdc42 and FgRac1 were involved in sexual development and pathogenesis while FgRho2 and FgRho4 were both involved in cell wall integrity, only FgRho4 showed a role in nuclear division and septum formation. The knocking down of TaRab7 enhanced the susceptibility of wheat Suwon 11 to an avirulent race CYR23, which implies that TaRab7 plays an important role in the early stage of wheatstripe rust fungus interaction and in stress tolerance [28]. In our research the Unigene17170_Tc15_2 peaked at 12 hpi in THTT and peaked at 48-72 hpi in THTS when forming appressorium, it may be involved in the formation of hyphae.
qRT-PCR revealed more detail differential expressed patterns at early stage between THTT

and THTS
The results of the qRT-PCR were almost same as the patterns in the RNA-seq at 144 hpi. However, there are more differential expressed patterns at early stage. Twelve tested genes have obvious peak expression before 144 hpi (Fig 3).
Expression profiles of 9 genes found in THTT precedes that in THTS. CL3900.Contig1 peaked at 6 hpi both in THTT and THTS. CL3499.Contig1 peaked at 6 hpi in THTT and 36 hpi in THTS. CL6956.Contig1 peaked at 6 hpi in THTT, which is much higher than that of THTS. Expression profiles of 4 genes found in THTS precedes that in THTT. CL622.Contig3 expressed at 6 hpi, and the quantity of CL622.Contig3 was higher than that in THTS. Unigene22186 had two peaks-12 hpi and 144 hpi in THTS while one peak at 144 hpi in THTT. Unigene11935 peaked at 12 hpi and 36 hpi in THTS with more than 100 times of the quantity in THTT. Unigene14763 had almost same expression pattern in two races.
Unigene18727 was found with conserved domain of cytochrome P450 52A6. CYP53 family members play a key role in fungal colonization of plant material by detoxification of anti-fungal compounds released by plants or generated during plant material degradation. Moreover, CYP53 family members play a role in the generation of a secondary metabolite, veratryl alcohol, which is crucial in the degradation of the plant cell wall component, lignin [29]. The expression of Unigene18727_Tc15_2 peaked in THTT at 12 hpi and gradually increases and peaks at 72 hpi in THTS, Further, it hardly expresses at later stage. Combining with the function of cytochrome P450, we assume that the gene helps infect host. And this gene expressed earlier and stranger than that in THTS.
CL2376.Contig1 is a member of ATP binding cassette transporter C family. The ABC transporter belongs to a large ancient protein family Energy produced by ATP binding and hydrolysis is used for the ABC transporter participation in the substrate transport process, such as the RNA translation and transmembrane process that is required for DNA repair. Although the mechanism is unclear, it plays an important role in enhancing the ability of the pathogen to resist adverse external environment.
CL2376.Contig1_Tc15_2 had two peaks, 24 hpi and 48 hpi in THTS, and had only one peak at 12hpi in THTT, and the expressed quantity was significant higher than that in THTS, so it possibly plays part in pathogenic process. Chitinase can degrade most fungal cell walls, prevent or interrupt fungal infection, colonization and expansion in plants. Chitin deacetylase (CAD) belongs to chitinase. The enhanced CGA activity in the process of the formation of appressorium in Uromyces and appressorium formation failure in the polycarbonate (PC) artificial culture medium caused by CAD deletion mutants in Magnaporthe oryzae both indicate that CAD is involved in the formation of appressorium [30]. Furthermore, CAD has a promoting effect on the formation of fungal fruiting body.
CAD isolated from the basidiomycetes of Flammulina velutipes specifically express during the fruiting stage . [31] CL3499.Contig2_Tc15_2 was predicated to encode chitin deacetylase. This gene expressed up regulatedly and peaked at 6 hpi in THTT and peaked at 72 hpi in THTS, expressed more earlier than that in THTS. Therefore, we speculated that the CL3499.Contig2_Tc15_2 gene may play a key role in the pathogenic infection of leaf rusts.
Effector Unigene22186 is with superoxide dismutase (SOD) activity. Jeong et al. [32] found that MnSOD role in organisms facing exogenous oxidative ROS and especially superoxide overproduction.
The role of MnSOD was deduced from the SOD2 increased expression after exposure to menadione in S. pombe. Furthermore, SOD2 mutants were more sensitive than wild type to menadione, and plumbagin, a menadione derivative. SOD2 plays an important role in the virulence of both C. neoformans var. gattii and C. neoformans var. grubii, depending on the route of infection [33,34].
Since defense against ROS is determinant for pathogenicity, MnSOD evolutionary history and pathophysiological roles in invasive or allergic mycoses agree with the hypothesis that pathogenicity emerged multiple times within fungi emerged multiple times within fungi. MnSOD was used for taxonomic and evolutionary data [35]. This gene expressed in THTT and up regulated from 6hpi and peaked at 144 hpi, however it was expressed in THTS had two peaks at 12hpi and 144 hpi respectively and the expressed quantity was higher than that in THTT. We guess that Unigene22186 is a disease-related gene.
There are many factors that lead to the difference of the two strains. In order to explore the reasons for the toxicity difference of these two strains, we will focus on these specific expression genes in future experiments.

Effectors secreted by THTT and THTS when interacted with Thatcher have different expression patterns
Functional domains were obtained by Interpro Scan to make clear the biological functions of the 54 candidate secreted proteins (Table S6) Candidate effector Unigene15605 was up regulated in THTS. It was found with Kre9/Knh1 family related to cell wall organisation in fungi. Gilbert et al. [36] have identified that the H99 deletion mutants kre5Δ and kre6Δskn1Δ contained lessβ-1,6-glucan, which is the component of cell wall, grew slowly with an aberrant morphology. Moreover, these two mutants resulted in alterations in cell wall chitosan and the exopolysaccharide capsule, a primary cryptococcal virulence determinant.
25-66AA of Unigene11683 is similar with concanavalin A-like lectin family. Sequence analysis were conducted to show that the N-terminal regions of BEACH proteins shares similarities with concanavalin A (ConA)-like lectin superfamily while members of the BEACH family are generally defined as trafficking regulatory of vesicle, a transport carrier in the process of secretory protein [37].
However, 21-69AA of the gene is aligned with glycoside hydrolase family 7. Van et al. [38] constructed knock-down (KD) mutants of cellulases in order to clarify that cellulases, which can hydrolyze crystalline cellulose to permeate the host epidermis belong to glycosyl hydrolase (GH) families 6 and 7. The results showed that transcript levels of the target genes and cellulase activity were considerably reduced in the KD mutants, and the KD mutants resulted in fewer lesions, less penetration, and infection of fewer cells compared with the parent strain.
Fifty-four candidates differentially expressed effectors and expression profiles of 4 candidate effectors were obtained and two candidate effectors were identified in this study. These results lay the foundation for the study of the pathogenic differences in Pt at the molecular level and reveal the interaction mechanisms. Furthermore, this study is mainly based on bioinformatics, macroscopic screening and qRT-PCR technology, and the specific differentially expressed genes should be further investigated. Other members of our experimental group have performed studies using gene silencing techniques and tobacco inoculation of the pathogenic genes to verify their function.

Conclusion
The results obtained in this study might have a solid foundation for the future studies on clarifying the mechanism of pathogenicity differences of Chinese isolates and pathogenic mechanism of Pt. The pathogenesis of wheat leaf rust is very complex although the pathotypes were very similar. The virulence results in significant changes at the molecular level in only a small number of genes according to the virulence formula. Therefore, exploring the pathogenic mechanism is a long-term, arduous and crucial task. Consistently exploring the wheat leaf rust effectors, understanding their function in pathogenesis, and locating and targeting their receptors in plants are indispensable. We will focus on several effectors that differentially expresssed in THTT and THTS to explore the avirulent genes. Therefore, this study lays the foundation for the understanding of the molecular mechanisms of wheat resistance, the control of wheat rust, the resistance to persistent diseases and the enrichment of the pathogenic mechanism of specific parasitic fungi.  Table 1 The      Student's t-test. P value < 0.05, *; P value < 0.01, **. n=4.

Figure 4
The candidate effector proteins that can inhibit the cell death triggered by Bax. A and B represent Unigene23118 and Unigene17565, respectively. C is the schematic of injection.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.