Plant response to infection
The only differences between strain THTS and strain THTT were the “S” and “T” in the corresponding pathotypes, which meant that 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).
DEGs based on RNA-seq
The total numbers of reads obtained from the databases for THTS and THTT were 7,206,771 and 7,040,346, respectively. The gene expression levels were compared between THTT as the treatment group and THTS as the control group. A total of 21,172 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. In addition, 2,784 significant DEGs (FDR ≤0.001 and a difference ≥2 times), of which 1,708 were upregulated and 1,076 were downregulated (Fig. 1A), were obtained. In total, 45 and 26 DEGs were expressed specifically in THTT and THTS, respectively. (Fig. 1B). Few of them were annotated including zinc metalloprotease, phospholipase, RNA polymerase, cell wall protein, actin cytoskeleton-regulatory complex protein and glycoprotein glucosyltransferase among the upregulated genes (Table S3), while phosphate acyltransferase, dihydroneopterin aldolase and adenosine triphosphatase (ATPase) among the downregulated genes (Table S4).
In THTT vs THTS, 2,784 DEGs were enriched for terms in the three GO functional categories: biological processes, cellular components and molecular functions (Fig. 2A and Table 2). A total of 875 genes were annotated to the biological process category, 539 genes were annotated to the cellular component category, and 971 genes were annotated to the molecular function category.
In the KEGG pathway analysis, a Q value less than or equal to 0.05 was taken to indicate that a DEG was significantly enriched. In total, 1,562 genes from among 2,784 DEGs were annotated to 153 KEGG pathways. The pathways with the highest enrichment confidence, 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, and MAPK signaling pathway - yeast pathways, are shown in Fig. 2B.
Candidate effectors among the 2,784 DEGs
One hundred fifty-nine proteins containing signal peptides were predicted to exist among the 2,784 proteins encoded by the DEGs using SignalP 4.1. A total of 118 proteins were found by TMHMM 2.0 to lack transmembrane domains. Then, EffectorP 2.0 [22] was used for further clarification, and 54 candidate effectors were ultimately screened. There were 21 upregulated candidate effector genes and 33 downregulated candidate effector genes in THTT compared with THTS. Structural analysis was conducted by a bioinformatics method (Table S5). The protein sizes ranged from 59 to 237 aa. In all, 41 candidate effector sequences contained more than three cysteines. The predicted 54 effector proteins were analyzed using the Nr database. The results showed that only Unigene22186 had a functional annotation, which was copper/zinc superoxide dismutase (SOD), and the rest were all hypothetical proteins. Domain analysis was performed using Pfam software, and glutaredoxin, DPBB_1, thaumatin, Calc_CGRP_IAPP, MF_alpha_N, TIG, Mtd_N, Omp_AT, Glyco_hydro_7, RNA_pol_Rpb2_3, GPI-anchored, CAP, OSTbeta, BAF250_C and Cys domains were found. Perl software was used to search for known motifs, including RxLR in oomycetes, [Y/F/W]xC in powdery mildew and G[I/F/Y][A/L/S/T]R in flax rust. The motif search revealed that 2 effector proteins contained the RxLR motif and that 20 effector proteins contained the [Y/F/W] xC motif; however, no proteins with the G[I/F/Y][A/L/S/T]R motif were found.
qRT-PCR analysis of 12 DEGs
The RNA-seq results were similar to the qRT-PCR results at 144 hpi. The gene expression trends were related to the formation of structures after rust infection. By 12 hpi, the germ tube, appressorium and substomatal vesicle were clearly formed. Primary hyphae and haustorial mother cells were formed by 18-24 hpi, and secondary hyphae and haustoria were formed by 36 hpi. By 48 hpi, small mycelial masses were formed [23]. To confirm that the trends in gene expression were consistent with the result of RNA-seq and to elucidate the expression characteristics at different time points, we selected important DEGs for real-time fluorescence quantitative analysis (Table 2).
Three genes enriched in different KEGG pathways were chosen for analysis. CL622.Contig1 (Fig. 3A) expression peaked at 288 hpi in THTT and at 144 hpi in THTS. CL3900.Contig1 (Fig. 3B) was significantly upregulated at 6 hpi in THTT and THTS (to a significantly greater extent in THTT than in THTS); after 6 hpi, the expression decreased. This gene was annotated as an ATPase. Based on the expression properties, the spores were in the stages of germination and bud tube formation at 6 hpi, so we speculate that CL3900.Contig1 is an energy carrier that participates in these processes. Unigene18070 (Fig. 3C) was annotated to be involved in helicase activity and DNA binding. According to the qRT-PCR results, this gene reached the highest expression level in THTS at 24 hpi, while it reached the highest expression level in THTT at 12 hpi, and the expression level in THTS was higher than that in THTT. Huang et al. [23] observed that Pt forms appressoria beginning at 8 hpi. Additionally, haustorial mother cells and primary hyphae develop when Pt contacts an epidermal or a mesophyll cell [24]. We suggest that appressorium formation occurred earlier in THTT than in THTS.
According to RNA-seq, Unigene1676 and Unigene17170 were specifically expressed in THTT at 144 hpi, while Unigene18727, CL3499.Contig2 and CL2376.Contig1 were specifically expressed in THTS at 144 hpi. However, according to qRT-PCR, these genes were not specifically expressed at time points other than 144 hpi. The peaks in expression of both Unigene1676 (Fig. 3D) and Unigene17170 (Fig. 3E) occurred at 12 hpi in THTT, earlier than those in THTS. Unigene17170 expression had two peaks at 48 hpi and 72 hpi in THTS. However, Unigene17170 expression peaked at 12 hpi in THTT. These genes were annotated as a ubiquitin (Ub) E3 ligase, a GTPase and an NADPH oxidase. Unigene18727 (Fig. 3F), CL3499.Contig2 (Fig. 3G), and CL2376.Contig1 (Fig. 3H) were expressed earlier in THTT than in THTS, and their expression peaked at 12 hpi in THTT. In contrast, in THTS, Unigene18727 expression peaked at 72 hpi, CL3499.Contig2 expression peaked at 36 hpi, and CL2376.Contig1 expression peaked at 24 hpi and 48 hpi. Unigene18727, CL3499.Contig2 and CL2376.Contig1 were annotated as cytochrome P450, ATP binding cassette transporter and chitinase, respectively.
The expression features of two candidate effectors, Unigene22186 and CL6956.Contig1, were analyzed. The expression trends of Unigene22186 (Fig. 3I), which exhibits SOD activity, were similar in both THTS and THTT, although Unigene22186 expression in THTS was higher than that in THTT. Two peaks were observed at 12 hpi and 144 hpi in THTS, and the highest peak appeared at 144 hpi. In THTT, the expression began to increase at 48 hpi and peaked at 144 hpi. The expression levels of CL6956.Contig1 (Fig. 3J), annotated to have hydrolase activity, were higher in THTT than in THTS. The expression of this gene peaked at 12 hpi in both pathotypes and then gradually decreased.
Two other differentially expressed effectors were also analyzed. The expression of Unigene11935 (Fig. 3K) was very low in THTT at every time point but peaked at 36 hpi in THTS. There was no expression of Unigene11683 (Fig. 3L) in THTT at 144 hpi. However, two peaks were detected at 12 hpi and 36 hpi. Unigene11683 expression peaked at 48 hpi in THTS. Therefore, Unigene11935, which was expressed only in THTS, was a specifically expressed DEG.
Unigene23118 and Unigene17565 could inhibit programmed cell death (PCD) induced by Bax in N. benthamiana
To test whether the effector proteins could inhibit PCD, the recombinant plasmids, Bax, and negative control-eGFP were injected into N. benthamiana. There was clearly no homologous recombination (HR) at the sites where the vector, Unigene23118 and Unigene17565 were injected alone. Necrosis occurred at the sites where Bax was injected alone and where Bax and the vector were injected together. However, there was no necrosis at the sites where Bax and either Unigene23118 (Fig. 4A) or Unigene17565 (Fig. 4B) were injected together, indicating that these genes inhibited Bax action, preventing PCD from being induced. Thus, we can come to a conclusion that Unigene23118 and Unigene17565 have effector functions.