RNAseq overview
We aimed to determine the changes in the N. benthamiana transcriptome that occur preceding cell death driven by non-host recognition of three (Zt9, Zt11, and Zt12) Z. tritici SSPs, and a secreted phytotoxic RNase (Zt6) in comparison to a green fluorescent protein (GFP) control. GFP, Zt6 and SSPs were transiently expressed in N. benthamiana leaves using agroexpression and samples were collected at 24- and 48-hours post-inoculation (hpi), i.e. prior to the HR becoming visible by eye. RNA was extracted from treated leaves to produce 30 RNA-seq libraries (five treatments x three biological replicates x two timepoints) for sequencing by paired-end sequencing on the Illumina Hiseq 2000 platform. The libraries contained 23.95-35.42 million raw reads. Subsequent quality filtering reduced the number of reads in each library by 37.4% to 51.15%. Of the remaining reads, 92.8% to 98.9% were successfully mapped to the reference N. benthamiana genome [10] (Table 1).
Principle component analysis (PCA) of the overall gene-expression profile showed that replicates of each treatment clustered tightly, as well as revealing minimal difference between the GFP and SSP treatments at 24hpi (Fig.1A). This pattern changed by 48hpi, with SSP treatments showing clear separation compared to the GFP control, although with minimal difference among themselves. In contrast, Zt6 expression induced a different gene-expression profile compared to both the GFP and SSP treatments at 24hpi. This was further exaggerated by the 48hpi timepoint (Fig.1A).
To examine these different responses, differential expression analysis was conducted using DESeq2 Bioconductor package [11], with FDR less than 0.05, to identify differentially expressed genes (DEGs) between the SSP and Zt6 treatments compared to the GFP control (Fig.1B). At 24hpi across all treatments, more DEGs were upregulated than downregulated. At 48hpi, the ratio of upregulated and downregulated DEGs was similar, although the total number of DEGs greatly increased. The increase in number of DEGs was most noticeable for the Zt9 treatment, which had a much smaller number of DEGs at 24hpi than the other treatments (Fig.1B).
The profound difference between the Zt6 and SSP treatments observed by PCA (Fig.1A), was reflected in the number of DEGs shared between treatments (Fig.2). For example, the number of upregulated genes common to all three SSP treatments at 48 hours (1294) was far higher than that of any other group of three treatments at 48hpi that contained Zt6 (732, 131, 127 DEGs respectively, Fig.2A). A similar pattern was observed for downregulated genes, although one particular Zt6-containing group of treatments (Zt11/Zt9/Zt6 vs GFP) at 48 hpi shared a much larger number of DEGs (1069) than other groupings, and almost as much as the SSP treatments group (Zt9/Zt11/Zt12 vs GFP, 1104 DEGs, Fig.1B). Noticeably, at 48 hpi the Zt11 treatment shared far more DEGs with the Zt6 treatment (482 and 808 upregulated and downregulated respectively) than either of the other two SSP treatments, but at 24hpi was much more similar to the Zt12 treatment (Fig.1).
Gene Ontology enrichment analysis
In order to examine whether differentially expressed genes (DEGs) were involved in specific developmental processes, we performed a GO-enrichment analysis of up- and down-regulated genes for the Zt6 (Tables 2-3) and SSP treatments (Tables 4-5). For brevity, only the top 20 GO terms (i.e. those with the largest number of DEGs averaged across treatments) enriched in the SSP treatments are listed here.
Most significantly enriched GO terms were present across several treatments, and most were enriched only at the 48hpi time point. Among upregulated DEGs, only four GO terms were enriched in more than one treatment at 24hpi, whilst none were enriched among downregulated DEGs. Out of all significantly enriched GO terms (70 in downregulated DEGs and 113 in upregulated DEGs) 37 were enriched only at the 24hpi time point of only one treatment. Noticeably, 20 of these were enriched only in the Zt12 treatment, while eight, seven, and two were exclusive to the Zt11, Zt9, and Zt6 treatments respectively. Only one of the eight GO terms exclusive to Zt11 (sulfate reduction) was enriched in the downregulated DEGs.
The contrasting transcriptional response to SSP and Zt6 treatments revealed by PCA (Fig.1A) was also observed in terms of enriched GO categories. Among the upregulated DEGs, only seven GO terms were enriched in all treatments at 48hpi (structural constituent of ribosome, translation initiation factor activity, intracellular, ribosome, translation, metabolic process, and ribosome biogenesis) (Tables 2,4). Most other GO terms were enriched in either one or all SSP treatments, or they were exclusive of the Zt6 treatment. There were only two terms (hydrolase activity, and catalytic activity) that were common to Zt6 and at least one other treatment. This pattern was similar in downregulated DEGs, where only nine GO terms were enriched across all treatments (carbohydrate metabolic process, extrinsic component of membrane, fructose 1,6-bisphosphate 1-phosphatase activity, membrane, photosynthesis, light harvesting, photosystem I reaction centre, photosystem II, and photosystem II oxygen evolving complex) (Tables 3,5) also in this case at 48hpi. Strikingly, these nine common GO terms included five (of a total six) photosynthesis-related GO terms. Of those remaining, only five were shared between the Zt6 treatment and at least one SSP treatment. Furthermore, whilst four photosynthesis related GO terms were also significantly enriched among upregulated DEGs, three were exclusive to the Zt11 treatment at 24hpi, and one was significantly enriched in both Zt11 and Zt12 treatments, also only at 24hpi.
The response to Zt6 expression was characterised by an overrepresentation of up-regulated genes involved in protein phosphorylation and kinase activity (Table 2). This included the GO terms transferase activity (of phosphorous-containing groups), protein phosphorylation, protein kinase activity, and protein serine/threonine kinase activity. By contrast, these signaling-related GO terms were enriched amongst downregulated genes for the SSP treatments (Table 5).
GO categories overrepresented among upregulated genes in the SSP treatments were those related to cellular protein catabolic processes. This included threonine type endopeptidase activity, endopeptidase activity, proteasome complex, and proteolysis involved in cellular catabolic process, as well as categories involved in other cellular protein metabolic processes such as protein folding (Table 4).
Over-represented among genes downregulated in the SSP treatments were those involved in diverse processes such as response to hormone and signal transduction, binding and activity of microtubules, lipid and fatty acid metabolism, DNA replication, and protein complex genes (e.g. MCM complex and kinesin complex) (Table 5). Closer investigation showed that downregulation of histone protein genes was the sole driver of enrichment of the GO term nucleosome. GTPase activity and GTP-catabolic process were the only GO terms enriched among both up and downregulated genes of the SSP treatments.
Differential expression of immune receptor-like genes
We previously demonstrated the requirement for the cell-surface co-receptors NbBAK1 and NbSOBIR1 for full induction of cell death by the SSP effector group [8]. These co-receptors are not required for Zt6-induced cell death [9]. It is therefore likely that cell death induced by the SSP treatments is a form of immune receptor-mediated programmed cell death. In contrast, Zt6-induced cell death is likely to be more similar to necrosis. Expression of receptors is often upregulated in response to the presence of their ligand [12, 13, 14]. We therefore assessed whether there were differential expression patterns of cell surface and cytoplasmic receptor gene families commonly associated with pathogen-associated molecular pattern (PAMP) and/ or effector recognition.
There were clear differences in transcriptional profiles across receptor families induced by Zt6 and SSP treatments. At 24hpi, transcriptional changes among wall-associated kinases (WAKs), receptor like kinases (RLKs), leucine-rich repeat receptor like-kinases (LRR-RLKs), and lectin-receptor kinases (LecRLKs) showed a clear bias toward upregulation in the Zt6 treatment (Fig.3). For WAKs and RLKs this bias was reversed in the SSP treatments, and for LRR-RLKs there was a smaller number of genes upregulated in the SSP treatments in comparison to Zt6. For LecRLKs, a similar bias towards upregulation existed between the Zt6, Zt11, and Zt12 treatments, although this was less prominent for Zt11 and Zt12. At 48hpi, LRR-RLKs in the SSP treatments showed a bias towards downregulation, this was less prominent in the Zt6 treatment. Interestingly, at 48hpi cysteine-rich receptor-like kinases (CRKs) showed a strong bias toward downregulation across all treatments. Genes annotated as nucleotide-binding site leucine-rich repeats (NLRs) showed bias toward upregulation at 48hpi across all treatments, though this was a comparatively small number of genes relative to the total N. benthamiana NLR gene complement.
Given the differences between the transcriptomes of the Zt6 and SSP treatments, we next focused on specific genes within these five receptor gene families that showed marked difference in expression change between the Zt6 and SSP treatments. We defined this as a log2 fold change of an absolute value of at least 1.0 in one or all of the SSP treatments which was either not present or reversed in the Zt6 treatment, or a log2 fold change of an absolute value of at least 1.0 in the Zt6 treatment which was either not present or reversed across all the SSP treatments. For all six receptor gene families, few expression changes were apparent at 24hpi in comparison to the GFP control. However, at 48hpi a clear difference in expression pattern is visible between the Zt6 and SSP treatments. Zt6 treatment specifically induced expression of 11 out of a total 149 LRR-RLKs, 4 out of a total 65 NLRs, 7 out of a total 107 WAKs/RLKs and 7 out of a total 54 LecRLKs in the current N. benthamiana genome annotation (Fig.4). Zt6 also induced expression of 34 out of a total 795 RLPKs in the current N. benthamiana genome annotation (Fig.5). Only two LRR-RLKs, two LecRLKs and 7 RLPKs were downregulated in response to Zt6 treatment (Fig.4). Interestingly, no CRKs showed any transcriptional response to Zt6 (Fig.4). Very few receptor genes were induced by the SSP treatments, although notable exceptions included one NLR gene induced by all three SSP treatments at 48hpi (Fig.4), and two RLPK genes relatively strongly induced only by the Zt12 treatment at 48hpi (Fig.5). Most receptor family genes that were transcriptionally responsive to the SSP treatments were downregulated. These included a notable over-abundance of WAKs and RLKs that were strongly downregulated in response only to the Zt11 treatment at 48hpi. Together, this suggests that there is reprogramming of receptor gene expression following exposure to SSPs, and that this differs between individual SSPs.
Differential expression of transcription factors, senescence- and programmed cell death-associated genes
There are differences in how developmental, pathogen-associated, and stress-induced PCD is executed. These distinct but partially overlapping pathways share some common components, although no pathway is fully characterised [15, 16, 17]. Our GO enrichment analysis showed that Zt6 and SSP treatments appear to be inducing different types of transcriptional reprogramming, and that Zt6-induced cell death is distinct from ordered receptor-mediated PCD [15]. We therefore investigated whether gene families commonly associated with transcriptional reprogramming related to senescence and PCD displayed differential expression patterns between the Zt6 and SSP effector treatments. We first evaluated the expression of the NAC, WRKY, TCP and MYB transcription factor (TF) families (Fig.6), which are widely reported as regulators of leaf senescence [18, 19, 20, 21]. At 24hpi there were few differential expression changes across any of the treatments. However, by 48hpi there were clear differences in expression patterns between Zt6 and SSP treatments. Overall, Zt6 treatment induced expression of 29 out of a total 241 NACs, 20 out of a total 151 WRKYs, 5 out of a total 61 TCPs and 18 out of a total 242 MYBs in the current N. benthamiana genome annotation. Only four TFs (all NAC family) were downregulated in response to Zt6 at this timepoint. In contrast, few TFs of any family were induced by the SSP treatments. Indeed, for WRKY, MYB and TCP families the majority of differentially expressed TFs were downregulated in response to effector treatment. There was minimal overlap between the individual genes transcriptionally induced by Zt6, and those transcriptionally repressed by the SSP treatments.
Several other non-TF gene families have previously been linked to senescence and cell death in plants. Metacaspases, zinc-finger domain containing proteins and HR-inducing proteins have all been shown to regulate PCD or HR [22, 23]. The senescence-associated genes SAG1 (Arabidopsis thaliana), SAG102 (Medicago truncatula), phytoalexin–deficient 4 (PAD4/SAG101) and harpin-induced gene 1 (HIN1), are markers of senescence and HR respectively [24, 25, 26]. In addition, STAY-GREEN (SGR) genes operate downstream of NAC TF regulation of senescence in the catabolism of chlorophyll [27]. We investigated expression of these genes across all treatments and found that expression of genes encoding metacaspases, SAG, HR-inducing and SGR proteins between Zt6 and SSP treatments was similar (Fig.S1). However, a zinc-finger containing protein encoding gene was strongly induced by Zt6 treatment at 24- and 48hpi but not by the SSP treatments.