
==== Front
Plant Commun
Plant Commun
Plant Communications
2590-3462
Elsevier

S2590-3462(24)00124-X
10.1016/j.xplc.2024.100882
100882
Research Article
Decomposition of dynamic transcriptomic responses during effector-triggered immunity reveals conserved responses in two distinct plant cell populations
Liu Xiaotong 1236
Igarashi Daisuke 14
Hillmer Rachel A. 17
Stoddard Thomas 1
Lu You 18
Tsuda Kenichi 15
Myers Chad L. 23
Katagiri Fumiaki katagiri@umn.edu
13∗
1 Department of Plant and Microbial Biology, University of Minnesota - Twin Cities, St Paul, MN 55108, USA
2 Department of Computer Science and Engineering, University of Minnesota - Twin Cities, Minneapolis, MN 55455, USA
3 Bioinformatics and Computational Biology Graduate Program, University of Minnesota - Twin Cities, Minneapolis, MN 55455, USA
4 Institute for Innovation, Ajinomoto Co., Inc., Kawasaki, Japan
5 State Key Laboratory of Agricultural Microbiology, Hubei Hongshan Laboratory, Hubei Key Lab of Plant Pathology, College of Plant Science and Technology, Huazhong Agricultural University, Wuhan 430070, China
∗ Corresponding author katagiri@umn.edu
6 Present address: 225 S 6th St, Minneapolis, MN 55402, USA

7 Present address: Biotalys, Inc., 2520 Meridian Parkway Suite 480, Durham, NC 27713, USA

8 Present address: Life Sciences Institute and Department of Cell and Developmental Biology, University of Michigan Medical School, Ann Arbor, MI 48109, USA

16 3 2024
12 8 2024
16 3 2024
5 8 10088215 8 2023
25 1 2024
13 3 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Rapid plant immune responses in the appropriate cells are needed for effective defense against pathogens. Although transcriptome analysis is often used to describe overall immune responses, collection of transcriptome data with sufficient resolution in both space and time is challenging. We reanalyzed public Arabidopsis time-course transcriptome data obtained after low-dose inoculation with a Pseudomonas syringae strain expressing the effector AvrRpt2, which induces effector-triggered immunity in Arabidopsis. Double-peak time-course patterns are prevalent among thousands of upregulated genes. We implemented a multi-compartment modeling approach to decompose the double-peak pattern into two single-peak patterns for each gene. The decomposed peaks reveal an “echoing” pattern: the peak times of the first and second peaks correlate well across most upregulated genes. We demonstrated that the two peaks likely represent responses of two distinct cell populations that respond either cell autonomously or indirectly to AvrRpt2. Thus, the peak decomposition has extracted spatial information from the time-course data. The echoing pattern also indicates a conserved transcriptome response with different initiation times between the two cell populations despite different elicitor types. A gene set highly overlapping with the conserved gene set is also upregulated with similar kinetics during pattern-triggered immunity. Activation of a WRKY network via different entry-point WRKYs can explain the similar but not identical transcriptome responses elicited by different elicitor types. We discuss potential benefits of the properties of the WRKY activation network as an immune signaling network in light of pressure from rapidly evolving pathogens.

This study reports a reanalysis of Arabidopsis transcriptome data collected during an effector-triggered immunity response, which reveals double-peak time-course patterns prevalent among upregulated genes. These double peaks, derived from single peaks in two cell populations, indicate a conserved immune response that is likely mediated by activation of the WRKY transcription factor network. Involvement of this network suggests potential adaptability in mounting resilient immune responses against rapidly evolving pathogens.

Key words

Arabidopsis
Pseudomonas syringae
plant immune response
multi-compartment model
WRKY transcription factors
resilient immune signaling network
Published: March 16, 2024
==== Body
pmcIntroduction

Two modes of inducible immunity have been well characterized in plants: pattern-triggered immunity (PTI) and effector-triggered immunity (ETI) (Jones and Dangl, 2006; Dodds and Rathjen, 2010). PTI signaling is initiated when molecular patterns in the apoplastic space are recognized by pattern-recognition receptors (PRRs) on the cell surface (Albert et al., 2020). Such molecular patterns include microbe-associated molecular patterns (MAMPs) and damage-associated molecular patterns (DAMPs). For example, a 22-amino-acid peptide of the bacterial flagellin, flg22, is a MAMP recognized by an Arabidopsis PRR, flagellin-sensitive 2 (FLS2) (Gómez-Gómez and Boller, 2000). Similarly, Arabidopsis DAMPs, plant elicitor peptides (PEPs), are recognized by the functionally redundant Arabidopsis PRRs PEP receptor 1 (PEPR1) and PEPR2 (Krol et al., 2010; Yamaguchi et al., 2010). Pathogens that are well adapted to plant hosts deliver effectors into the plant cell that interfere with PTI signaling (Khan et al., 2018). ETI signaling is initiated when some of the pathogen effectors or their effects are recognized by cognate intracellular resistance (R) proteins (van der Hoorn and Kamoun, 2008). For example, the effects of AvrRpt2 and AvrRpm1 effectors from the bacterial pathogen Pseudomonas syringae on the plant protein RIN4 are recognized by the Arabidopsis R proteins resistance to Pseudomonas syringae 2 (RPS2) and resistance to Pseudomonas syringae maculicola 1 (RPM1), respectively (Mackey et al., 2002, 2003; Axtell and Staskawicz, 2003).

Plant immunity is expensive, as it requires large amounts of energy and resources (Bolton, 2009). Furthermore, plant immunity includes responses that are harmful to host cells as well as pathogens. For example, during ETI, plant cells that autonomously recognize a pathogen effector undergo programmed cell death, called a hypersensitive response (HR) (Coll et al., 2011). The cost of immunity is likely the reason why plants have evolved inducible, rather than constitutive, immunity (Heidel et al., 2004). It is important that plants rapidly induce immunity at the right place when it is needed. Thus, tracking when and where particular immune responses are induced with respect to pathogen infection is crucial for understanding plant immunity.

Transcriptome analysis is a standard, relatively economical approach for measurements of many different responses. For example, Mine et al. collected time-course transcriptome data during ETI in six Arabidopsis genotypes inoculated with P. syringae pv. tomato DC3000 (Pto) strains carrying an empty vector or expressing AvrRpt2 or AvrRpm1 effectors (Mine et al., 2018). This study revealed the importance of immune signaling sectors (Tsuda et al., 2009) in rapid transcriptome responses during ETI. There are other studies in which transcriptome analysis was performed at the cellular level of spatial resolution during plant immunity (e.g., Chandran et al., 2010; Zhu et al., 2023). However, transcriptome analysis with sufficiently high resolution in both time and space is technically and economically challenging, particularly with plant tissues.

Time-course transcriptomes are typically analyzed assuming independence across time points using profile-based approaches such as clustering or statistical models such as ANOVA (Eisen et al., 1998; Ben-Dor et al., 1999; Cui and Churchill, 2003). Although these approaches are useful for classification of genes into different expression-pattern groups and detection of differentially expressed genes, respectively, the temporal continuity of the data is not fully appreciated. Mathematical models that can capture temporal continuity information, such as regression and spline models, have also been used (Ramsay and Silverman, 2005). However, these models are not built upon specific mechanistic hypotheses and generally fail to provide mechanistic interpretations.

A multi-compartment model (MCM) is a type of mathematical model used to describe some quantity (signals, in our study) transmitted among compartments of a system (Godfrey, 1983). It has been applied to many fields of biology, including pharmacokinetics and epidemiology (Kermack and McKendrick, 1927; Rescigno, 1960; Prabakaran et al., 2021). Each compartment has an input and an output and a rule defining the output. A simple implementation of a single compartment has a first-order self-decay as the rule (Figure 1A). The output of a compartment (compartment Y in Figure 1A) is specified by the input (output of compartment X), the input amplification ratio (aY), and the self-decay rate (kY). An MCM formed by a network of such simple compartments is a tractable linear ordinary differential equation (ODE) system, and the mechanistic parameters of each compartment, the input amplification ratio and self-decay rate, are readily interpretable.Figure 1 Two populations of plant cells, spatially defined relative to the locations of bacterial cells, can explain the double-peak time courses of mRNA levels.

(A) A two-compartment model. Each of the compartments, X and Y, has a parameter that determines the first-order self-decay rate (kx and ky). The input level to compartment Y is proportional to the input amplification rate parameter ay. The output of X has a first-order decay, whereas the output of Y has a single-peak time course (see the time courses shown above the X and Y boxes). The model is defined by ordinary differential equations (ODEs) (1) and (2) for each compartment. The analytical solutions for the output of the compartments are generally in the form of linear combinations of first-order decays in MCM. Solutions for the two-compartment model under the initial conditions of X=1 and Y=0 at t=0 are shown.

(B) The MCM structure used in this study. A series of signaling compartments (S(1) to S(6), shown as an example) generate single-peak time-course patterns with delayed peak times (shown above each compartment) when S(n)(0)=0,n=1,2,…,6. The signaling compartment series is common among all the genes we modeled, and the parameter values for the signaling compartments were predetermined. For each gene, the output from a signaling compartment (S(1) as an example in the figure) was used as the input to the first-peak (ACP) response compartment (RA,g), and the output from another signaling compartment (S(6) as an example in the figure) was used as the input to the second-peak (NACP) response compartment (RN,g). A common gene-specific decay rate, kg, was assumed for the ACP and NACP response compartments of a gene g. The gene-specific input amplification ratio parameters of the response compartments are aA,g. and aN,g. Thus, for each gene, the selection of the two signaling compartments and the parameters kg, aA,g, and aN,g were fitted to the mRNA-ratio time-course data corresponding to log2(Pto AvrRpt2/mock).

(C) A schematic diagram of the two plant cell populations. The plant cells in the close vicinity of bacterial cells receive the effector AvrRpt2 from the bacterial cells, and the ETI response is induced cell autonomously (ACP; blue cells). Most cells in the ACP undergo HR cell death by 9 hpi. The cells surrounding the ACP respond to AvrRpt2 in an indirect manner (NACP; orange cells). NACP cells do not die.

(D) A schematic diagram explaining how a single-peak response in each of the ACP and NACP (blue and orange traces, respectively) can be observed as a double-peak response (red trace) in the log(mRNA-level ratio) time course for gene g. The mRNA molecules from two cell populations are mixed when the tissue is macerated during the RNA preparation procedure.

(E) We decomposed the double-peak time-course pattern into two single-peak patterns (first and second peaks) using MCM. Each peak can be characterized by the peak level, peak width, and peak time.

In this study, we reanalyzed the transcriptome dataset generated by Mine et al. (2018). We noticed that double-peak time-course patterns were prevalent among the genes upregulated in Arabidopsis leaves upon inoculation with a relatively low dose of Pto AvrRpt2. We implemented an MCM to decompose the double-peak time course into two single-peak time courses for each gene. This double-peak decomposition revealed an “echoing” pattern between the two single peaks in approximately 1400 upregulated genes: when the first peak for a gene was early, its second peak was also early; and when the first peak for a gene was late, its second peak was also late. We demonstrated that the first- and second-peak responses likely represent the responses of two distinct plant cell populations: the cells responding cell autonomously to AvrRpt2 (autonomous cell population [ACP]) and those surrounding the ACP and responding indirectly (non-autonomous cell population [NACP]). The echoing pattern of transcriptome responses also showed that there is a well-conserved transcriptome response between the two cell populations, whereas their eliciting signal-molecule types are different. We found that genes highly overlapping with the conserved gene set were also upregulated during PTI elicited by flg22, with dynamics similar to the ACP response during ETI. The echoing mRNA response pattern between the two cell populations appeared to be regulated mainly by WRKY transcription factors. We propose WRKY network activation through different entry-point WRKYs as a mechanism to explain the similar but not identical transcriptome responses across the ACP and NACP during ETI and PTI. Architectural features of the proposed signaling network centered on the WRKY activation network suggest network properties potentially beneficial for an immune signaling network that is under strong pressure from rapidly evolving pathogens.

Results

Double-peak time-course patterns are prevalent in the mRNA levels of genes upregulated by Pto AvrRpt2

From the dataset (Mine et al., 2018), we selected 3039 genes significantly upregulated at one or more time points after 3 h post inoculation (hpi) with Pto AvrRpt2 compared with the mock treatment (water). We focused on the time range after 3 hpi because major responses of most genes were confined to this time range (Supplemental Figure 1A). We further selected 2435 high-precision upregulated genes (Supplemental Figure 1B and 1C). A preliminary visual inspection using a spline model of the mRNA-level log-ratio of Pto AvrRpt2 over mock for each gene suggested widespread double-peak time-course patterns (Supplemental Figure 2A).

We hypothesized that widespread double-peak responses are the superimposed responses of two distinct cell populations (Lu and Tsuda, 2021), each of which has a single-peak response (Figure 1C–1E). In fact, the experimental system with a low inoculation dose of Pto AvrRpt2 (optical density 600 [OD600] = 0.001; Mine et al., 2018) contains at least two plant cell populations with distinct responses (e.g., Hatsugai et al., 2016). One population of plant cells, the ACP, receives and cell autonomously responds to the ETI-eliciting effector AvrRpt2. Most ACP cells undergo HR cell death by 9 hpi (Hatsugai et al., 2016). Another population of plant cells, the NACP, surround the ACP cells and do not undergo HR cell death. Based on the cell death timing, we reasoned that any responses with peak times after 9 hpi must be attributed to NACP cells. Thus, a better-defined form of the two-cell population hypothesis is that responses with peak times before 9 hpi are largely attributable to ACP cells.

The MCM used in the study

To test the hypothesis of two cell populations, we designed an MCM to decompose the double-peak mRNA-level log-ratio pattern for each gene into two single-peak patterns. A chain of compartments generates a single-peak time-course pattern at each compartment with a delayed peak time along the chain (signaling compartments labeled as S(i); i=1,2,…,6 in Figure 1B). This signaling compartment chain may be considered to represent molecular components in a conceptualized signal transduction pathway. Because it models behaviors of two cell populations, this conceptualized signal transduction pathway includes a step for intercellular signaling between ACP and NACP cells; the steps before and after this intercellular signaling step represent signal transduction within ACP and NACP cells, respectively. For practicality of model fitting, the parameters of the signaling compartments, including their numbers, were predetermined.

For regulation of a particular gene g, we modeled as if the outputs from some upstream and downstream signaling compartments (S(1) and S(6) in Figure 1B as an example) were used in the ACP or NACP, respectively. The mRNA-level log-ratios of gene g in the ACP and NACP were modeled as the outputs of the response compartments RA,g and RN,g, respectively. When the leaf tissue was macerated for RNA preparation, mRNA molecules from the ACP and NACP were mixed, and the mixed mRNA-level log-ratio was modeled as the double-peak pattern, RA,g+RN,g. The input amplification parameters aA,g and aN,g of the response compartments depend on regulation of gene g in the ACP and NACP and the ratio of the ACP and NACP in the tissue. RA,g and RN,g represent two decomposed single-peak time courses for the gene (Figure 1B and 1E).

The MCM explains the double-peak transcriptome response well

To avoid bias potentially introduced by using predetermined signal-compartment parameters, consistency in the model fit between two sets of the parameter values, sets 1 and 2 (Text S1; the fitted parameter values are in Supplemental Table 1) was used as the criterion for a good model. The fitted values for the two peaks with sets 1 and 2 across 2435 selected upregulated genes were generally well correlated (Supplemental Figure 3). We selected 1889 genes with Pearson correlation ≥0.9 as the genes stably modeled by the MCM, and the MCMs with set 1 for these genes were used in subsequent work.

The MCM compared with the data over time for each gene is given in Supplemental Figure 4. The median across the genes of the Pearson correlations between the mean estimates and the MCM-modeled values across the time points was 0.93 (Supplemental Figure 5). The heatmaps in Supplemental Figure 6A show the discrepancies across the genes between the mean estimates of the data and the MCM-modeled values at each time point. The discrepancies were small and showed no obvious pattern, except for the 191 genes shown at the bottom of the figure. These 191 genes had very early first and second peaks, and their third peaks were not captured by the model. We concluded that our MCM explains the double-peak transcriptome response well.

The second-peak response echoes the first-peak response

Using the outputs of the response compartments, RA,g and RN,g, the double-peak pattern was decomposed into two single-peak patterns for each gene (Supplemental Figures 4 and 6B). We characterized the single-peak patterns by their peak times and peak levels (Figure 1E). We classified the genes into three groups, 214 first-peak-specific, 1366 echoing, and 227 second-peak-specific genes, according to the ratio of the first-peak level over the second-peak level (1807 modeled-and-classified upregulated genes; Supplemental Figure 7 and Supplemental Table 2), sorted the genes in each group according to their peak times, and visualized them using heatmaps (Figure 2A). Among the 1366 echoing genes, it was strikingly evident that the first-peak time-course patterns were well replicated in the second-peak time-course patterns: when the first peak of a gene was relatively early, its second peak tended to be early; when the first peak was late, its second peak tended to be late. We called this replicated pattern an echoing of transcriptome responses (hence the gene-group name). The second-peak responses were time stretched as well as delayed compared with the first-peak responses in the echoing genes (note the size differences in time intervals in Figure 2A). The mean ratio of the second-peak time over the first-peak time, which indicates the mean time-delay rate, was 2.8. The mean ratio of the second-peak width over the first-peak width (width at 75% peak level; Figure 1E), which indicates the time-stretching ratio, was 2.2 (Supplemental Figure 8). Because these values are similar, the second-peak response is approximately time-scaled from the first-peak response with a scaling factor of about 2.5.Figure 2 The second-peak response echoes the first-peak response across about 1400 genes.

(A) Heatmaps of the first-peak and second-peak responses modeled by MCM (left and right panels, respectively) for 1807 modeled-and-classified upregulated genes. The genes were classified into three groups, first-peak (ACP)-specific, echoing, and second-peak (NACP)-specific genes, according to the peak-level-ratio value between the first and second peaks of the gene. The number of genes in each group is also shown. Within each group, the genes were ordered according to the modeled peak time values of the first peaks of the genes for first-peak-specific and echoing genes and to those of the second peaks for second-peak-specific genes. The values used in the heatmaps are the log2 mRNA-level ratio of Pto AvrRpt2-inoculated over mock at the indicated hpi. Among the 1889 modeled upregulated genes, 82 genes were unclassified and are not included in this figure.

(B) A heatmap of the flg22 response mean estimates (Hillmer et al., 2017). The gene groups and the gene orders are the same as in (A). The log2 mRNA-level ratio of flg22-infiltrated Arabidopsis Col-0 over the flg22-infiltrated Arabidopsis fls2 mutant was normalized to the vector size of 1 for each gene across the time points. These normalized log2FC values were used in the heatmaps. Although the original data included a 1-hpi time point, those data were removed from the figure because such early responses largely represent non-immune-specific, general stress responses (Bjornson et al., 2021).

We interrogated our hypothesis that the first-peak responses, whose peak times were mostly earlier than 9 hpi, are largely attributable to the ACP. First, we compared the transcriptome response after Pto AvrRpt2 inoculation to that after Pto inoculation in the same dataset (Mine et al., 2018). With Pto AvrRpt2 inoculation, ACP cells received the ETI-eliciting effector AvrRpt2, in addition to immune-suppressing effectors, and underwent the autonomous ETI response, including the HR. With Pto inoculation, the cells spatially equivalent to ACP cells received only immune-suppressing effectors, and their immune response was thus predicted to be strongly suppressed during the comparable time range. Consistent with this prediction, in the time range for the first-peak response in the Pto-AvrRpt2-inoculated samples (4–9 hpi), very little response was observed in the Pto-inoculated samples (Figure 1D in Mine et al., 2018; Supplemental Figure 9). This observation suggests that the transcriptome responses in the early time range are mainly responses of cells spatially equivalent to the ACP.

Second, we studied RNA sequencing (RNA-seq) data from Arabidopsis Col-0 plants in which AvrRpt2 protein was transgenically and conditionally expressed in plant cells (in-planta-expressed AvrRpt2 data) (Hillmer et al., 2023). Because AvrRpt2 protein was expressed upon induction in virtually all cells, resulting in macroscopic HR cell death (Tsuda et al., 2012; Hatsugai et al., 2017), the transcriptome response before cell death in the in-planta-expressed AvrRpt2 data should mainly represent responses in the ACP. When the peak times and peak levels of the first- and second-peak responses to Pto AvrRpt2 and in-planta-expressed AvrRpt2 were compared, the first-peak response was more similar than the second-peak response to the response to in-planta-expressed AvrRpt2 (Figure 3A–3D). Next, we sorted 1807 modeled-and-classified upregulated genes according to the ratio of the first-peak level over the second-peak level, and subsets of the sorted genes were obtained using a sliding window of 150 genes. When the numbers of overlapping genes between each of the gene subsets and the 1972 genes upregulated by in-planta-expressed AvrRpt2 were visualized, the overall trend was clear. The higher the ratio of the first-peak level over the second-peak level, the greater the number of overlapping genes (Figure 3E). Because genes with higher peak ratios were more first-peak specific/preferential, this trend also indicated that the first-peak response is more similar to the response to in planta-expressed AvrRpt2. In summary, two lines of evidence support the hypothesis that the first-peak response is largely from the ACP.Figure 3 The first-peak response is more similar than the second-peak response to the ACP response induced by in-planta-expressed AvrRpt2.

(A–D) Peak times (A and C) and peak levels (B and D) of the first-peak (A and B) and second-peak (C and D) responses were compared to peak times (A and C) and peak levels (B and D) of the response to in-planta-expressed AvrRpt2 across the 1092 (A and B) and 967 (C and D) overlapping genes. Because late-peak-time estimates for the in-planta-expressed AvrRpt2 response are not precise, genes with in-planta-expressed AvrRpt2 peak times later than 7.5 hpi are aggregated into a group, and a boxplot of the Pto AvrRpt2-inoculated peak times for the group is shown (>7.5 hpi; A and C). A linear regression line (black line) and Pearson correlation coefficient (PCC) are shown in each panel. In (A) and (C), the linear regression and the correlation coefficients were calculated for the genes excluding the >7.5-hpi group after the peak time values were offset subtracted and log transformed, as seen in the scaling of the axes.

(E) The level of overlap in the modeled upregulated genes between the Pto AvrRpt2-inoculated samples and the in-planta-expressed AvrRpt2 samples along the ratio of peak levels between the first and second peaks in the Pto AvrRpt2-inoculated samples is shown using a sliding window method (window size, 150 genes). The dashed horizontal line and the gray shading indicate the mean and its 95% confidence interval (Fisher’s exact test, two-sided), respectively, for all 1807 modeled-and-classified upregulated genes. Pink vertical dashed lines show the boundary first-peak level/second-peak-level ratio values of four gene groups: ACP-specific, echoing 1, echoing 2, and NACP-specific genes. Echoing genes shown in Figure 2A were divided into echoing 1 and echoing 2 genes at the first-peak-level/second-peak-level ratio for the upper boundary of the 95% confidence interval.

Two-peak decomposition sorted genes into functionally distinct groups

If gene groups defined by their two-peak decomposition characteristics are significantly enriched in genes with different functions, we could conclude that our two-peak decomposition captures characteristics of actual biological mechanisms. We tested this idea using Gene Ontology (GO)-term enrichment analysis. As classifiers, we used the first- and second-peak level ratio and the peak times. There were four peak-level-ratio categories (ACP-specific, echoing 1 [ACP preferential], echoing 2 [NACP preferential], and NACP-specific) (Figure 3E) and two peak-time categories, early and late, for a total of 4 × 2 = 8 groups. Note that the early and late peak times are determined relatively within each peak-level-ratio category and are not directly comparable across peak-level-ratio categories.

The results of biological process GO-term enrichment analysis are summarized in Figure 4 (details in Supplemental Table 2). General defense-related terms (“defense response to …”) were enriched in ACP-preferential and early groups, but this trend was weaker in the echoing 2 late group, and no defense-related terms were significantly enriched in NACP-specific groups (Supplemental Table 2). This observation may suggest that the ETI transcriptome responses specific to the NACP are quite different from well-characterized immune-related transcriptome responses.Figure 4 Enrichment of biological process GO terms specific to gene groups defined by double-peak characteristics.

The values in the gene-group columns show fold enrichment compared with the genome mean proportions. The higher the fold enrichment value, the darker the orange color. The sizes of the eight gene groups are shown in the second row. The “whole.genome” column shows the number of genes in the genome for each GO term.

HR-related terms were specifically enriched in the echoing 1 early group, consistent with the fact that the ACP undergoes the HR. (We postulate that genes with HR-related terms are more likely to be marker genes associated with processes leading to the HR rather than HR executors.) These terms were not significantly enriched in the ACP-specific early group. The trend that HR-related genes tended also to be expressed, albeit relatively weakly, in the NACP suggests that the HR program could be activated in the NACP as well as the ACP by default but that an HR-suppressing program could also be activated in the NACP to limit spread of the HR (Jabs et al., 1996). HR-associated accumulation of reactive oxygen species (ROS) (Su et al., 2018) appears to stress both the ACP and NACP, as the “response to oxidative stress” term was enriched in gene groups for both cell populations. Hypoxia response terms were strongly enriched in the ACP-specific early group. These observations suggest that O2 is converted to ROS mainly in the HR-undergoing ACP, and O2 is thus limiting in the ACP, inducing the hypoxic response. ACP-generated ROS probably diffuse and stress both the ACP and NACP. The “sulfate assimilation” term was enriched in the echoing 2 late group, suggesting that reducing compounds such as glutathione accumulate in the NACP, perhaps as part of the HR-suppression program (Király et al., 2012).

Terms related to the water-deprivation response, including response to abscisic acid, were enriched in the echoing 1 early and late groups and the echoing 2 and NACP-specific late groups. Limiting water accessibility in the apoplastic space is an effective ETI response against apoplastic pathogens such as Pto (Xin et al., 2016). Plants may activate a water-deprivation response to draw water from the apoplast into the cytosol, thereby limiting apoplastic water availability, rapidly in the ACP and later in the NACP.

It is common for the mRNA levels of photosynthesis-related genes to be downregulated during immune responses (Bilgin et al., 2010) (as in this dataset; Supplemental Table 3), which would decrease the synthesis of new photosynthetic apparatus. Enrichment of chlorophyll-catabolic genes in the NACP-specific early group suggested active downregulation of the photosynthetic apparatus specifically in the NACP. Terms related to fatty-acid catabolism were enriched in the echoing 2 late group, suggesting mobilization of an alternative energy source in the NACP while photosynthesis was rapidly downregulated. In summary, gene grouping based on two-peak decomposition characteristics revealed group-specific/preferential gene ontology (GO)-term enrichment, demonstrating the biological relevance of our two-peak decomposition approach.

Echoing genes with binding sites for WRKYs, NACs, HSFs, and CAMTA1 are associated with certain peak times

To gain insight into the regulatory mechanisms underlying the echoing transcriptome response, we performed transcription factor (TF)-binding site-enrichment analysis of 1807 modeled-and-classified upregulated genes. For 38 highly significantly enriched TF-binding sites, we visualized the distributions of over- or under-representation of genes with each TF-binding site along the peak times for first- and second-peak responses (Figure 5A). Most of the TFs for the binding sites in the figure belong to three large TF families, the NAC (No Apical Meristem/Arabidopsis thaliana Transcription Activator Factor 1/2/Cup-shaped Cotyledon 2; Mathew and Agarwal, 2018), WRKY (tryptophane-arginine-lysine-tyrosine; Birkenbihl et al., 2018; Chen et al., 2019), and HSF (heat shock factor; Guo et al., 2016) families. The time-dependent representation patterns of the WRKY and HSF families were similar, whereas there appeared to be two different patterns within the NAC family (the heatmaps and the dendrogram on the left in Figure 5A). Whereas members of both NAC groups (NAC1s and NAC2s, blue and cyan TF names) were overrepresented around 8 hpi, only NAC2s were overrepresented around 6.5 and 14 hpi and under-represented around 20 hpi. WRKYs (orange TF names) showed a common pattern in the first- and second-peak responses: strong overrepresentation around 4.5 and 10 hpi, followed by overrepresentation around 6.3 and 18.5 hpi, and then by strong under-representation around 8 and 20 hpi. Thus, the representation patterns of genes with WRKY-binding sites showed an echoing pattern. HSFs (red TF names) were overrepresented around 5 and 16.5 hpi, but the number of sliding windows with overrepresentation in the first-peak response was higher, suggesting greater overrepresentation in ACP-specific genes. Indeed, the proportions of the genes with binding sites for any of the HSF members were 26%, 16%, and 11% for ACP-specific, echoing, and NACP-specific genes, respectively. This higher representation of genes with HSF-binding sites among ACP-specific genes is consistent with the heat response GO-term enrichment in the ACP-specific early gene group (Figure 4).Figure 5 Genes with WRKY-binding sites are abundant among 1366 echoing genes and show an echoing distribution pattern across the first- and second-peak responses.

(A) Enrichment of genes with the indicated TF-binding sites along the peak time of the first- and second-peak responses. The significance in a −log10 scale and whether the genes with the TF-binding sites are over-represented or under-represented are shown in the heatmaps. The color scale bar is shown at the bottom. Members of the NAC1, NAC2, WRKY, and HSF families are indicated by blue, cyan, orange, and red font colors. CAMTA1 and dof24 are singletons and are shown in black. Structurally, AT3G42860 and AT3G09735 do not belong to these TF families, but they were included in the WRKY and HSF families, respectively, owing to their similar distribution patterns along the peak times. The horizontal bar plot on the right shows the proportion of genes with the indicated TF-binding sites in each of the first- and second-peak-response gene sets (214 ACP-specific plus 1366 echoing genes and 1366 echoing genes plus 227 NACP-specific genes, respectively).

(B) Heatmaps similar to (A) with consolidated NAC1s, HSFs, and WRKYs and CAMTA1 for 1366 echoing genes only. The horizontal bar plot on the right shows the proportion of echoing genes with the TF-binding sites.

We consolidated each of the TF families into one for better visualization of singletons, calmodulin-binding transcription activator (CAMTA)1 and DNA binding with one finger (dof) 24: for one consolidated TF family, genes with binding sites for any of the family members were counted. The consolidated TF-binding-site heatmap generally conserved the features of the TF-family patterns described above (Supplemental Figure 10).

To further explore the echoing transcriptome response, we selected NAC1s, HSFs, WRKYs, and CAMTA1, which showed significant temporal variation in representation (Figure 5B). Genes with the CAMTA1-binding site had an early echoing overrepresentation around 4.5 and 10 hpi. CAMTAs are transcriptional activators for general stress responses, likely responding to increased intracellular Ca2+ concentration (Doherty et al., 2009; Benn et al., 2014; Bjornson et al., 2021). An increase in intracellular Ca2+ concentration also occurs in the ACP early during ETI (Grant et al., 2000), which explains the overrepresentation of genes with the CAMTA1-binding site around 4.5 hpi. Their overrepresentation around 10 hpi strongly suggests that elicitation of responses in the NACP is also associated with increased intracellular Ca2+ concentration. Because the proportion of echoing genes with CAMTA1-binding sites was only 9% and because its echoing pattern was limited to early times, CAMTA1 cannot explain the response behaviors of most echoing genes. In this light, WRKYs were remarkable: 68% of echoing genes had WRKY-binding sites (the significance of the WRKY-binding site enrichment by Fisher’s exact test, two-sided, was p = 9.2 × 10−84), and the echoing representation patterns of the genes cover most of the time range of interest (Figure 5B). We concluded that it is highly likely that WRKYs are central to the echoing transcriptome response during ETI.

Diverse WRKY genes are transcriptionally upregulated in the time ranges of interest in both cell populations

We examined transcriptional regulation of 35 WRKY genes among the 3039 genes upregulated after Pto AvrRpt2 inoculation. Their temporal mRNA levels are shown in a heatmap (Figure 6A), and the profiles of binding sites for the WRKY proteins in the dataset (columns) (O'Malley et al., 2016) for the 35 WRKY genes (rows) are also shown (Figure 6C). Different subgroups of WRKYs (Eulgem et al., 2000) are indicated by different colors of WRKY gene or protein names. We speculate that binding specificities of the WRKY proteins are more similar within a subgroup than across subgroups. First, many WRKY genes were transcriptionally upregulated from 4 to 24 hpi, strongly suggesting that many WRKY TFs were activated in the time ranges of interest in both cell populations (Figure 6A). Second, there was no obvious correlation between the mRNA-level time-course pattern and certain WRKY subgroups (i.e., the colors of the WRKY gene names are well mixed in Figure 6C), suggesting that some members of most WRKY subgroups were upregulated during the time ranges for both cell populations. This means that any genes with at least one WRKY-binding site could be transcriptionally regulated by WRKYs at any point during the time range in either cell population. Third, 28 of the 35 upregulated WRKY genes had at least one WRKY-binding site (Figure 6C), suggesting that the majority of the 35 WRKY genes were themselves regulated by WRKY TFs. Although WRKY TFs can be transcriptional activators or repressors (Birkenbihl et al., 2018), these observations are consistent with the notion that WKRY TFs likely play a major role in mediating the echoing transcriptome responses.Figure 6 Diverse WRKY genes, most of which have WRKY-binding sites, are transcriptionally activated during the times of interest.

(A) Heatmap of the log mRNA-level ratio (Pto AvrRpt2 over mock) time courses for 35 WRKY genes after inoculation with Pto AvrRpt2.

(B) Heatmap of the log mRNA-level ratio (wild type over fls2) time courses for 35 WRKY genes after infiltration with flg22.

(C) Numbers of binding sites for diverse WRKY TFs (labeled at the bottom) in the 35 WRKY genes (labeled at the right). Different font colors for the names of WRKY genes (right) or TFs (bottom) represent different subgroups defined by Eulgem et al. (2000). The WRKY gene orders in the row are the same for (A) to (C). The color scale bars for log2FC in (A) and (B) and for the number of binding sites per gene in (C) are shown at the right of (C).

(D) A schematic diagram of the WRKY-network activation hypothesis. The signals initiated by distinct elicitors on the left, AvrRpt2 (an ETI elicitor), PEP (a DAMP), and flg22 (a MAMP), are fed into a transcriptional activation network of WRKY TFs (WRKY Network in the middle) through different entry-point WRKYs. Different colors of WRKY letters represent WRKY subgroups with somewhat different binding specificities. The elicitors and WRKYs that are activated by the elicitors are indicated by the same colors of enclosing boxes (green, pink, and orange for AvrRpt2, PEP, and flg22). Note that, although different elicitors do not activate the same set of WRKYs, every elicitor activates members of most WRKY subgroups. Non-WRKY TFs (e.g., TF-A, shown at the bottom of the WRKY network) can be part of the network if their genes have WRKY-binding sites, and some WRKY genes have binding sites for the non-WRKY TFs. WRKY-controlled genes on the right have somewhat different WRKY-binding sites (W-boxes). Different W-box colors represent the binding sites for the WRKY subgroups of the same colors. Thus, as long as some members of most subgroups are activated, very similar sets of WRKY-controlled genes are upregulated, which explains the similar but not identical transcriptome responses across the different elicitors. WRKY-controlled non-WRKY TFs (e.g., TF-B at the bottom right) can expand the WRKY-controlled gene set.

Application of the discovery to interpretations of a single-nucleus RNA-seq dataset

After the bioRxiv preprint of the current study (Liu et al., 2022) was posted, a bioRxiv preprint of an Arabidopsis single-nucleus RNA-seq (snRNA-seq) study with biological samples related to the current study was posted by Nobori et al. (2023). We decided to analyze and interpret their snRNA-seq data using our new discovery of the transcriptomic responses of two distinct cell populations during ETI.

Analysis of snRNA-seq data is challenging because the total read counts per nucleus are typically limited (1077 reads per nucleus was the median of the Nobori et al. data subset we used). Therefore, it is not realistic to interpret snRNA-seq data at a “single genes × single nuclei” resolution. Furthermore, comparisons of snRNA-seq data to bulk RNA-seq data (the latter is dominated by reads from cytosolic poly(A)+ RNA), are complicated by the fact that nuclear and cytosolic poly(A)+ RNA profiles are not the same (Zaghlool et al., 2021). We began our analysis with “sets of many genes × single nuclei.” Based on our discoveries, we divided the genes in the Nobori et al. dataset into six non-overlapping gene sets: ACP early-up gene set; ACP middle-up gene set; ACP late-up gene set (the early, middle, and late genes were determined by the maximum expression time among 4, 6, and 9 hpi, respectively, in the first-peak response of the ACP-specific and echoing genes in Figure 5A), NACP-specific up gene set (Figure 5A), ETI Down gene set (Supplemental Table 3), and non-responsive gene set (348, 1008, 221, 226, 4813, and 26 084 genes, respectively). Note that the majorities of the ACP early-, middle-, and late-up genes were echoing genes and were also upregulated in the NACP at later time points. The ratio of all read counts for each of the first five responsive gene sets to those of the non-responsive gene set was calculated for each nucleus as the between-nuclei normalized gene-set expression value. The subset of snRNA-seq data we used from Nobori et al. consisted of two bacterial strains (Pto and Pto AvrRpt2) at three time points (4, 6, and 9 hpi), in addition to a mock inoculation treatment (9 hpi), with no replication. The bacterial inoculation dose (OD600 = 0.001) was the same as that used for the bulk RNA-seq dataset of Mine et al. (Mine et al., 2018). Thus, we expected to observe the responses of the two cell populations, ACP and NACP, in the Nobori et al. dataset. Additional details of the methods and results are described in Text S2.

We analyzed the distribution of normalized gene-set expression levels across the nuclei in each snRNA-seq sample for each responsive gene set. The plots in Figure 7A–7E are the quantile functions of the log2-transformed normalized gene-set expression levels in decreasing order. Any differences in the quantile function curves indicate that some fraction of the nuclei have different normalized gene-set expression levels between the snRNA-seq samples: a shift of part of the curve for a treated sample upward or downward from the mock curve indicates that some fraction of nuclei have the gene set expressed higher or lower in the sample than in the mock, respectively. With no replicates in the dataset, it was not possible to estimate the confidence intervals of the mock curves. However, because the curves for the Pto 4 and 6 hpi snRNA-seq samples, in which no strong transcriptome responses were expected (Mine et al., 2018), were very close to the mock curves for all five gene sets, we assumed that the confidence intervals of the mock curves were relatively narrow.Figure 7 The snRNA-seq data can be interpreted in the context of the two cell populations in the ETI transcriptome response.

(A–E) Quantile functions of the normalized gene-set expression levels of the snRNA-seq samples (the color and line type codes are shown on the right of E) in decreasing order for the indicated gene sets. (A) ACP early-up genes; (B) ACP middle-up genes; (C) ACP late-up genes; (D) NACP-specific up genes; (E) ETI down genes.

(F–H) The maximum cosine correlations across the genes in the color-coded gene sets (the color codes are in I) between the log-transformed normalized gene expression levels in the nuclear subset defined by a sliding window along the percentile of the nuclei (the x axis) in the indicated snRNA-seq samples and the MCM-modeled log-transformed expression-level time course for either the ACP or the NACP response (indicated).

(I–K) The MCM-model times for the maximum cosine correlations shown in (F) to (H), respectively. Dashed vertical lines represent the upper boundaries of the nuclei discussed in the text. Data from Nobori et al. (2023).

First, it was obvious that the NACP response was induced in almost all nuclei in the Pto AvrRpt2, 9 hpi snRNA-seq sample, as the entire curves for the ACP early-, middle-, and late-up gene sets were clearly above the mock curves (Figure 7A–7C) and the entire curve for the ETI down genes was clearly below the mock curve (Figure 7E). Note that, by this time, most of the ACP cells were dead, and nuclei from ACP cells should not have been represented in this snRNA-seq sample. We estimated the time correspondence in the NACP response between the Pto AvrRpt2 9 hpi snRNA-seq sample and our MCM-modeled NACP response based on bulk RNA-seq data. For this purpose, the profiles across many genes × subsets of many nuclei were analyzed. The nuclei of the snRNA-seq sample were assigned to subsets using a sliding window with 8% width (a sliding step of 1%) along the nuclear percentile (x axis) determined for the normalized gene-set expression levels of each of the ACP early-, middle-, and late-up gene sets (Figure 7A–7C). Within each subset of nuclei, the normalized expression level was calculated for each gene of the gene set. The profile of log2-transformed normalized gene expression level across the gene-set genes was compared to the profile of bulk-RNA expression-level predictions for the NACP response by the MCMs for the corresponding genes along the time course (NACP model time; from 3.4 to 18 hpi with a 0.2-h interval). Figure 7F shows the maximum cosine correlation for each sliding window across the model time for each gene set, and Figure 7G shows the NACP model times for the maximum correlations. The nuclear percentile values in the figures are the middle values of the sliding windows. The maximum correlation values and the model times for the maximum correlations were relatively stable across the nuclei. The average NACP model time across the nuclei was 13.4 hpi; i.e., the transcriptome response of most nuclei at 9 hpi of Pto AvrRpt2 inoculation approximately corresponded to the response of the cytosolic NACP transcriptome at 13.4 hpi (13.4 hpi/9 hpi ∼1.5-fold time scaled).

Second, we looked for some indication of the ACP response using earlier snRNA-seq samples with Pto AvrRpt2 inoculation. It appeared that the Pto AvrRpt2 4 hpi curve for the ACP middle-up genes was slightly above the mock curve toward their left ends (Figure 7B and its close-up view in Supplemental Figure 11B). We applied the sliding-window approach to the ACP middle-up gene set in the Pto AvrRpt2 4 hpi snRNA-seq sample with the ACP response MCMs. We found that the profiles up to the 15th nuclear percentile (window middle value) were relatively highly correlated, with the average of the corresponding ACP model times of 5.7 hpi (Figure 7H and 7I; 5.7 hpi/4 hpi ∼1.4-fold time scaled). This estimate of 15% of the cells as ACP cells is likely an underestimate, because it appears that nuclei from late-stage ACP cells are missing from the snRNA-seq data (see below). Considering this possible underestimation, the estimate of 15% ACP cells is reasonably close to 19%, which was estimated on the basis of microscopy observations of the HR dead cells following the same inoculation dose of Pto AvrRpt2 (Figure 2C; “wild type” of Hatsugai et al., 2016).

Third, the Pto AvrRpt2 6 hpi curve for the ACP early-up gene set was clearly above the mock curve toward their left ends (Figure 7A and its close-up view in Supplemental Figure 11A). Relative to the timing of the ACP response above, this must be the beginning of the NACP response. We applied the sliding-window approach to the ACP early-up gene set in the Pto AvrRpt2 6 hpi snRNA-seq sample with the NACP response MCMs. We found that the profiles up to the 17th nuclear percentile were relatively highly correlated, with the average of the corresponding NACP model times of 10.4 hpi (Figure 7J and 7K; 10.4 hpi/6 hpi ∼1.6-fold time scaled). The observation that only 17% of the nuclei initiated the NACP response at 6 hpi, whereas almost all nuclei showed the NACP response at 9 hpi (see above), suggests that the NACP is dynamically heterogenous; some cells begin the NACP response earlier than other cells. When we compared the nuclear responses across the gene sets, we observed such dynamic heterogeneity among the nuclei at 9 hpi as well (Figure TS2.2 in Text S2).

With the approximately 1.5-fold time scaling (the mean of three cases above) between the nuclear and cytosolic poly(A)+ RNA responses, we expected to see the ACP late-up gene set upregulated in the nuclei of the ACP cells in the Pto AvrRpt2 4- and/or 6-hpi snRNA-seq samples. However, we observed little upregulation of the ACP late-up gene set in either of the samples (Figure 7C and further analysis in Text S2). This lack of detection of the late-stage ACP response suggests that nuclei of late-stage ACP cells that highly expressed the ACP late-up genes were selectively removed from the snRNA-seq samples. Methods for nuclear isolation from plant green tissue typically include harsh treatments, such as high-concentration detergents or nuclear sorting, to remove chloroplasts, which are far more abundant than nuclei in crude preparations (Nobori et al. used a detergent treatment). It is conceivable that nuclei of late-stage ACP cells close to their death were more fragile and may have been easily damaged by such treatment. Damaged nuclei leak RNA and lose their representation in snRNA-seq data. It will be important to evaluate potential data biases introduced by nuclear isolation methods when applying snRNA-seq technologies to plant green tissues.

We used some of our discoveries to interpret the snRNA-seq data of Nobori et al. The discoveries used included (1) the existence of two cell populations, ACP and NACP, with distinct dynamic behaviors during the ETI response; (2) gene sets whose responses are defined dynamically and in a cell-population-specific manner; and (3) gene response MCM predictions for the ACP and NACP. We learned that the ACP cells are, at minimum, 15% of the leaf cells at a Pto AvrRpt2 inoculation dose of OD600 = 0.001 and that the NACP cells are dynamically heterogeneous. We also learned two characteristics of the snRNA-seq data. First, compared with the ETI response of cytosolic poly(A)+ RNA, that of nuclear poly(A)+ RNA is approximately 1.5 times faster. Second, nuclear poly(A)+ RNA information from late-stage ACP cells appears to be lost in the snRNA-seq data.

Discussion

MCM is a powerful modeling framework for temporal responses

MCM is a general mathematical framework that is highly applicable to a wide range of systems. The structure of the MCM is a network of compartments. Thus, MCM can be applied to any system that can be described by a directed graph. Any quantity flowing through the network can be modeled, e.g., concentrations of a drug in different parts of the body in pharmacokinetics (compartments represent the body parts) (Rescigno, 1960), numbers of people in different infectious stages in epidemiology (compartments represent the stages) (Kermack and McKendrick, 1927; Prabakaran et al., 2021), or intensities of signals in different parts of a regulatory network in our study (compartments represent the regulatory parts). Each compartment is defined by the input amplification ratio and self-decay rate parameters, and the behavior of each compartment is described by an ODE with the parameters. Typically, the ODEs describing an MCM are linear, and their analytical solutions can be expressed parametrically.

We designed an MCM structure with (1) a chain of signaling compartments, representing a conceptualized signaling pathway common for transcriptional regulation of all genes; and (2) two gene-specific response compartments, each of which draws a signal from a signaling compartment and represents the mRNA response of the gene in one of the two cell populations (Figure 1B). The outputs of the two response compartments were combined to represent the mixing of mRNA from two cell populations during mRNA sample preparation (Figure 1C and 1D). With parametrically described analytical solutions available, we fitted the analytical solutions to the data, as this was more computationally efficient than fitting the corresponding ODEs. This made it practical for us to fit the MCM to the RNA-seq read-count data, assuming a negative binomial residual distribution.

We predetermined the decay rate parameters of 11 signaling compartments in our MCM. This was because the data to constrain the MCM were available only for gene-specific outputs of the MCM and did not have sufficient power to constrain these parameters. We fit the outputs of which signaling compartments to use (a discrete variable), two input amplification ratio and one common decay rate parameters of two response compartments to the data for each gene. A consequence of predetermining the signaling compartment parameters was that the response compartment parameters were not mechanistically interpretable. This is the reason why we analyzed the kinetic characteristics we learned from the MCM (e.g., peak levels and times) instead. If more biological information about the signaling compartments were available, such as the output of each signaling compartment, we could build a more realistic model for the chain of signaling compartments, which would allow mechanistic interpretations of the response-compartment parameter values.

Elicitor(s) for the NACP transcriptome response likely originate from the responding ACP

Although recognition of the effector AvrRpt2 by the cognate R protein RPS2 through the guardee/decoy RIN4 initiates the signaling event inside the ACP cells (Axtell and Staskawicz, 2003; Mackey et al., 2003), it is not clear which elicitor molecules are recognized for initiation of signaling events in the NACP cells. There are two potential sources of elicitor(s) for transcriptome responses in the NACP, which are not mutually exclusive: signaling molecule(s) from the responding ACP and/or diffusible molecule(s) from the pathogen Pto AvrRpt2. Potential molecules derived from the bacterium are limited to diffusible molecules because the spread of bacterial cells from initial sites in the plant leaf is limited before 24 hpi (Zhu et al., 2023).

One candidate class of signaling molecules that could originate from the responding ACP cells are DAMPs (Lu and Tsuda, 2021). For example, PEP peptides are DAMPs whose precursors are encoded by PROPEP1 to 6 (Huffaker et al., 2006). The mRNA time-course pattern of PROPEP3 showed early echoing patterns (AT5G64905 in Supplemental Figure 4). Thus, both the ACP and NACP produced PEP very quickly. It is possible that the amount of PEP produced by ACP cells may not be sufficient to elicit all cells in the NACP. Thus, a wave of PEP may be radially propagated through the NACP from the ACP cells at the center. This wave of PEP propagation could explain the time-stretched nature of the echoing transcriptome response pattern in the NACP (Supplemental Figure 8B). The dynamic heterogeneity of NACP cells detected in the snRNA-seq dataset (Figure 7H and 7K and Text S2) was also consistent with a wave of PEP propagation.

DAMPs are a class of molecules that can induce PTI (Yamaguchi et al., 2010; Liu et al., 2013). Another major class of molecules that can induce PTI are MAMPs such as flg22 (Gómez-Gómez and Boller, 2000). We therefore compared the echoing transcriptome response patterns of the two cell populations upon Pto AvrRpt2 inoculation with the flg22-induced transcriptome response (Hillmer et al., 2017) (Figure 2A and 2B). Remarkably, the time-course patterns of the 1366 echoing genes identified in the Pto AvrRpt2 response were similar in flg22-induced responses (Supplemental Figure 12), although the flg22-induced responses were generally weaker and variable (for flg22 responses, the mRNA-level ratio values were normalized for each gene). The similarity in the time-course patterns of NACP and flg22 responses further supports the notion that the elicitor(s) for the Pto AvrRpt2 responses in the NACP may be DAMPs. Note that the actual peak times of the flg22 response, in which all cells received flg22 at 0 hpi, was more similar to the peak times of the ACP response than the NACP response, consistent with the PEP propagation wave hypothesis for the time-stretched response in the NACP. Furthermore, the similarity in the time-course patterns among all ACP, NACP, and flg22 responses implies that major portions of the transcriptome responses in ETI and PTI share a common regulatory mechanism.

Activation of the WRKY network may be a common signaling mechanism in the two distinct cell populations

What molecular mechanism could induce transcriptome responses with similar relative timing in two cell populations during ETI and in cells during PTI? We found that WRKY TFs may be important for such transcriptome responses in the ACP and NACP during ETI (Figure 5). Eighty percent of the WRKY genes upregulated upon Pto AvrRpt2 inoculation contained WRKY-binding sites (Figure 6C), suggesting that once some WRKY TFs are activated transcriptionally or post-transcriptionally, their encoded WRKY TF proteins bind and transcriptionally activate other WRKY genes. Many diverse WRKY genes form a transcriptionally interconnected network, and a substantial part of the WRKY network could eventually be transcriptionally activated once one or a few entry-point WRKYs are activated. Although the sets of WRKYs activated through different entry-point WRKYs are unlikely to be identical, such activated WRKY sets contain members of many WRKY subgroups, which could result in similar downstream transcriptome responses (Figure 6D). Any other TF genes with WRKY-binding sites and the TFs that can bind to some WRKY genes could also be part of the WRKY network (e.g., TF-A in Figure 6D). In fact, many WRKY genes have binding sites for diverse TFs (Supplemental Figure 13). Such TF-A and TF-B (see below) TFs could be responsible for upregulation of the approximately 30% of echoing genes that lack WRKY-binding sites.

Consistent with this WRKY-network activation hypothesis, we observed that at any time point from 4 to 24 hpi of Pto AvrRpt2, WRKY genes representing most WRKY subgroups were transcriptionally upregulated (Figure 6A). Activation of large portions of the WRKY network during this time range in the two cell populations suggests that any gene with some type of WRKY-binding site could potentially be upregulated in either of the cell populations. Because 68% of the 1366 echoing genes had WRKY-binding sites, this potential activation of the WRKY network by different elicitors could explain the similar gene sets for the echoing responses. Members of most WRKY subgroups were transcriptionally upregulated between 4 and 24 hpi of flg22, whereas the specific WRKY genes upregulated in ETI and PTI were not the same (Figure 6A and 6B). Thus, a similar but not identical set of upregulated genes during PTI can also be explained by the WRKY-network activation hypothesis (Figure 6D).

How could the conserved echoing dynamics in the two cell populations be achieved under this WRKY-network activation hypothesis? There are two possible mechanisms, which are not mutually exclusive: indirect activation through other TFs and/or intrinsic diversity among response genes. First, many of the upregulated genes with WRKY-binding sites encoded TFs themselves. Genes that can be upregulated by such TFs would be upregulated with a delay (indirect activation through other TFs; e.g., TF-B in Figure 6D). Some NAC1s may be examples of TF-Bs and accumulate with a delay, which could explain the role of NAC1s in very late responses in both cell populations (Figure 5). Second, the time course of transcriptional activation of a particular gene could be affected by characteristics of the gene itself, which we call intrinsic diversity among response genes. For example, transcriptional activation of a particular gene with a WRKY-binding site could be modulated by other TFs and/or epigenetic mechanisms, and such modulating mechanisms are conserved between the ACP and NACP. Note that there is no difference between cells in the ACP and NACP before inoculation. Differences in the mRNA degradation rate would also result in different peak times. In short, indirect activation and intrinsic diversity among response genes could explain the conserved response timing across the echoing genes.

The structural organization of the proposed WRKY-network-based signaling mechanism can provide an advantage to plants in the evolutionary race against rapidly evolving pathogens. First, the WRKY network is highly interconnected, which would provide resilience to the network against targeting by pathogen effectors. Second, the resilience of the WRKY network leaves only small parts of the immune signaling mechanism, namely the pathways specific to each signaling event, vulnerable to targeting by pathogen effectors. Such limited, vulnerable parts of the signaling mechanisms could be protected by, for example, the use of parallel signaling pathways connecting to different entry-point WRKYs or R proteins guarding the pathway. Third, the massive and conserved immune transcriptome response could be relatively easily repurposed or reutilized for new signaling pathways. For example, recognition of a new pathogen-associated molecule could be rapidly connected to the massive immune transcriptome response through evolution. Also, if a specific pathway were compromised by a pathogen effector, the signal could be rapidly rerouted through a different entry-point WRKY.

It should be emphasized that, in addition to this massive transcriptome response conserved across different types of immunity, there were also specific transcriptome responses. We observed upregulation of ACP- and NACP-specific genes (Figure 2). There were also upregulated genes specific to flg22-elicited PTI (Supplemental Figure 14). It appears that the massive, conserved transcriptome response provides a resilient core response across different types of immunity, while different types of immunity have further adapted to their specific needs.

In this study, we implemented an MCM to decompose the double-peak pattern, which was prevalent among the mRNA-level time courses of ETI-upregulated genes, into two single-peak patterns. We demonstrated that the decomposed first and second single-peak responses largely represent responses of two cell populations, the ACP and NACP, respectively. This temporal mRNA accumulation pattern was also relatively conserved during PTI. We propose that activation of the WRKY network is central to the generation of this massive and conserved immune transcriptome response across different types of immune responses. The properties of the proposed signaling network potentially provide benefits, as the immune signaling network is under strong pressure from rapidly evolving pathogens.

Materials and methods

Detailed methods are provided in Text S1.

RNA-seq data analysis, general

The data analysis was performed in R (R_Core_Team, 2021) using custom codes. All R scripts are available from https://github.com/fumikatagiri/MCM_At_AvrRpt2_time_course.

RNA-seq dataset, mRNA-level value estimates, and gene selection

We used a subset of the RNA-seq read-count dataset from Mine et al. (2018) (GEO: GSE88798). The subset included data collected from the leaf tissues of Arabidopsis Col-0 at 1, 2, 3, 4, 6, 9, 12, 16, 20, and 24 hpi of mock (water) or Pto AvrRpt2 with three biological replicates. We used the entire dataset, which contains 366 libraries, for the pre-selection of 18 442 genes that had at least 39 read counts for the 15th highest among 366 read-count values for each gene.

We fit a negative binomial generalized linear model (GLM-NB) to the read-count data, while normalizing the 90th percentile read count of each library (Hillmer et al., 2017). All mRNA levels and level ratios in this study are expressed as log2-transformed values.

We selected 3039 significantly upregulated genes by comparing the mean values of Pto AvrRpt2 samples to those of mock samples after 3 hpi (Storey FDR = 0.05, log2fold change [FC] >1; Supplemental Figure 1A). We defined the signal-to-noise ratio as (maximum of log2FC across the time points)/(mean of the SEs for log2FC across the time points) and selected 2435 high-precision upregulated genes with a signal-to-noise ratio >6.5 (Supplemental Figure 1B and 1C). A spline model was fit to the log2FC data for each of 2435 genes (Supplemental Figure 2).

An MCM for a double-peak time-course

The output change rate of a single compartment is defined by ODE (2) in Figure 1A. Our MCM was composed of a chain of 12 signaling compartments, S(n)(t),n=0,1,2,…,11, which generated single-peak time courses with delayed peak times. Two gene-specific response compartments, RA,g and RN,g, linearly combined the outputs of two single-peak signaling compartments with different peak times to generate a double-peak time-course pattern for gene g (Figure 1B). Two sets, sets 1 and 2, of predetermined decay rates for the signaling compartments were used to generate the single-peak patterns with the predetermined peak times at signaling compartments S(n),n=1,3,5,7,9,11: set 1, ([4,6,9,12,16,20]×0.95) hpi and set 2, ([3.5,5,7.5,10.5,14,18]×0.95) hpi, which were relatively orthogonal.

For each gene, the outputs of two signaling compartments were used as inputs to the response compartments, and the best of 10 combinations (indicated by the n value of S(n): (1, 5), (1, 7), (1, 9), (1, 11), (3, 7), (3, 9), (3, 11), (5, 9), (5, 11), and (7, 11)) was selected through MCM fitting by considering the combination as a discrete variable of the model. We assigned a common decay rate parameter, kg, to both response compartments, RA,g and RN,g, for each gene (Figure 1B) to avoid overfitting. Thus, for gene g, four parameters, kg, two input amplification ratios of the response compartments, aA,g and aN,g (Figure 1B), and the signaling compartment combination, were fit in the MCM, given one of sets 1 or 2 for the signaling compartment parameter values.

The MCM was fitted to the log-transformed mRNA ratio (log2(Pto AvrRpt2/mock), which was equal to log2(Pto AvrRpt2) – log2(mock)). This was because mRNA species in Arabidopsis commonly follow a double-exponential decay time-course pattern (Sorenson et al., 2018). Consequently, the mRNA-level log-ratio of a gene should follow a first-order decay, which is assumed at each compartment.

Altered MCMs for genes that had a time-course pattern of a single peak with a shoulder

The modeled values were clearly higher than the data values at 6 hpi for 191 genes (the gene selection criteria are given in Text S1), suggesting that they had early second-peak responses and that their time courses resembled a single peak with a shoulder. To impose early double peaks in model fitting, only the read-count data up to 16 hpi were used in the altered MCM for these 191 genes (Supplemental Figure 15).

TF-binding site enrichment analysis

We used the Arabidopsis cistrome dataset for 349 TFs across >30 000 genes in the analysis (O'Malley et al., 2016). For each gene, the DNA sequence 1000-bp upstream to 200-bp downstream of the transcription start site was defined as the promoter region and searched for TF-binding sites. Fisher’s exact test (one-sided) was used to calculate the p value for enrichment of genes with each TF-binding site among ACP-specific, echoing, or NACP-specific genes compared with the proportions in the genome. This yielded 100 TFs with a Benjamini-Hochberg FDR-corrected p < 0.01.

We examined whether the proportion of genes with a binding site for each of the 100 significant TFs showed temporal changes. Subsets of (ACP-specific + echoing) or (echoing + NACP-specific) gene sets were obtained using a sliding window of 150 genes with a sliding step of 15 genes along the peak time of the first-peak or the second-peak response, respectively. For each gene subset, the p value was calculated using Fisher’s exact test (two-sided) compared to the proportion in the gene sets. We selected 38 TFs with p < 0.001 in at least one window subset for further study. Similar analysis was performed for consolidated TF families and echoing genes only.

Funding

This work was supported by grants from the 10.13039/100000001 National Science Foundation (grant nos. MCB-0918908 and MCB-1518058 to F.K. and C.L.M. and IOS-1645460 to F.K.), a grant from the United States Department of Agriculture-National Institute of Food and Agriculture to F.K. (grant no. 2020-67013-31187 ), and a grant from Ajinomoto Co., Inc. to F.K. We thank the Minnesota Supercomputing Institute for their computing resources. We thank Tatsuya Nobori for information on the gene symbols in his snRNA-seq data.

Author contributions

X.L., C.L.M., and F.K. conceived the project, developed the MCM approach, analyzed the data, and wrote the manuscript. D.I., R.A.H., T.S., Y.L., and K.T. generated the unpublished RNA-seq dataset (Hillmer et al., 2023), which provided important supporting evidence for this study.

Supplemental information

Document S1. Supplemental Texts 1 and 2 and Supplemental Figures 1–15

Supplemental Figure 2A. Spline models of 2435 high-precision upregulated genes

Supplemental Figure 4. MCM with decompositions and data points for 1889 modeled upregulated genes

Supplemental Figure 15. Altered models were fit for 191 genes

Supplemental Table 1. The fitted parameter values in MCMs with set 1 or set 2 of signaling-compartment decay parameters

Supplemental Table 2. List of ACP-specific, echoing 1, echoing 2, and NACP-specific genes (each divided into early and late gene groups) and results of their GO-term enrichment analysis

Supplemental Table 3. List of ETI-downregulated genes (4848 genes) and results of their biological process GO-term enrichment analysis

Document S2. Article plus supplemental information

Acknowledgments

No conflict of interest is declared.

Published by the Plant Communications Shanghai Editorial Office in association with Cell Press, an imprint of Elsevier Inc., on behalf of CSPB and CEMPS, CAS.

Supplemental information is available at Plant Communications Online.
==== Refs
References

Albert I. Hua C. Nürnberger T. Pruitt R.N. Zhang L. Surface Sensor Systems in Plant Immunity Plant Physiol. 182 2020 1582 1596 10.1104/pp.19.01299 31822506
Axtell M.J. Staskawicz B.J. Initiation of RPS2-specified disease resistance in Arabidopsis is coupled to the AvrRpt2-directed elimination of RIN4 Cell 112 2003 369 377 10.1016/s0092-8674(03)00036-9 12581526
Ben-Dor A. Shamir R. Yakhini Z. Clustering gene expression patterns J. Comput. Biol. 6 1999 281 297 10.1089/106652799318274 10582567
Benn G. Wang C.Q. Hicks D.R. Stein J. Guthrie C. Dehesh K. A key general stress response motif is regulated non-uniformly by CAMTA transcription factors Plant J. 80 2014 82 92 10.1111/tpj.12620 25039701
Bilgin D.D. Zavala J.A. Zhu J. Clough S.J. Ort D.R. DeLucia E.H. Biotic stress globally downregulates photosynthesis genes Plant Cell Environ. 33 2010 1597 1613 10.1111/j.1365-3040.2010.02167.x 20444224
Birkenbihl R.P. Kracher B. Ross A. Kramer K. Finkemeier I. Somssich I.E. Principles and characteristics of the Arabidopsis WRKY regulatory network during early MAMP-triggered immunity Plant J. 96 2018 487 502 10.1111/tpj.14043 30044528
Bjornson M. Pimprikar P. Nürnberger T. Zipfel C. The transcriptional landscape of Arabidopsis thaliana pattern-triggered immunity Nat. Plants 7 2021 579 586 10.1038/s41477-021-00874-5 33723429
Bolton M.D. Primary Metabolism and Plant Defense-Fuel for the Fire Mol. Plant Microbe Interact. 22 2009 487 497 10.1094/mpmi-22-5-0487 19348567
Chandran D. Inada N. Hather G. Kleindt C.K. Wildermuth M.C. Laser microdissection of Arabidopsis cells at the powdery mildew infection site reveals site-specific processes and regulators Proc. Natl. Acad. Sci. USA 107 2010 460 465 10.1073/pnas.0912492107 20018666
Chen X. Li C. Wang H. Guo Z. WRKY transcription factors: evolution, binding, and action Phytopathol. Res. 1 2019 13
Coll N.S. Epple P. Dangl J.L. Programmed cell death in the plant immune system Cell Death Differ. 18 2011 1247 1256 10.1038/cdd.2011.37 21475301
Cui X. Churchill G.A. Statistical tests for differential expression in cDNA microarray experiments Genome Biol. 4 2003 210 10.1186/gb-2003-4-4-210 12702200
Dodds P.N. Rathjen J.P. Plant immunity: towards an integrated view of plant-pathogen interactions Nat. Rev. Genet. 11 2010 539 548 10.1038/nrg2812 20585331
Doherty C.J. Van Buskirk H.A. Myers S.J. Thomashow M.F. Roles for Arabidopsis CAMTA transcription factors in cold-regulated gene expression and freezing tolerance Plant Cell 21 2009 972 984 10.1105/tpc.108.063958 19270186
Eisen M.B. Spellman P.T. Brown P.O. Botstein D. Cluster analysis and display of genome-wide expression patterns Proc. Natl. Acad. Sci. USA 95 1998 14863 14868 10.1073/pnas.95.25.14863 9843981
Eulgem T. Rushton P.J. Robatzek S. Somssich I.E. The WRKY superfamily of plant transcription factors Trends Plant Sci. 5 2000 199 206 10.1016/s1360-1385(00)01600-9 10785665
Godfrey K. Compartmental Models and Their Application 1983 Academic Press
Grant M. Brown I. Adams S. Knight M. Ainslie A. Mansfield J. The RPM1 plant disease resistance gene facilitates a rapid and sustained increase in cytosolic calcium that is necessary for the oxidative burst and hypersensitive cell death Plant J. 23 2000 441 450 10.1046/j.1365-313x.2000.00804.x 10972870
Guo M. Liu J.H. Ma X. Luo D.X. Gong Z.H. Lu M.H. The Plant Heat Stress Transcription Factors (HSFs): Structure, Regulation, and Function in Response to Abiotic Stresses Front. Plant Sci. 7 2016 114 10.3389/fpls.2016.00114 26904076
Gómez-Gómez L. Boller T. FLS2: an LRR receptor-like kinase involved in the perception of the bacterial elicitor flagellin in Arabidopsis Mol. Cell 5 2000 1003 1011 10.1016/s1097-2765(00)80265-8 10911994
Hatsugai N. Hillmer R. Yamaoka S. Hara-Nishimura I. Katagiri F. The μ Subunit of Arabidopsis Adaptor Protein-2 Is Involved in Effector-Triggered Immunity Mediated by Membrane-Localized Resistance Proteins Mol. Plant Microbe Interact. 29 2016 345 351 10.1094/MPMI-10-15-0228-R 26828402
Hatsugai N. Igarashi D. Mase K. Lu Y. Tsuda Y. Chakravarthy S. Wei H.L. Foley J.W. Collmer A. Glazebrook J. Katagiri F. A plant effector-triggered immunity signaling sector is inhibited by pattern-triggered immunity EMBO J. 36 2017 2758 2769 10.15252/embj.201796529 28811287
Heidel A.J. Clarke J.D. Antonovics J. Dong X. Fitness costs of mutations affecting the systemic acquired resistance pathway in Arabidopsis thaliana Genetics 168 2004 2197 2206 10.1534/genetics.104.032193 15611186
Hillmer R.A. Igarashi D. Stoddard T. Lu Y. Liu X. Tsuda K. Katagiri F. The Kinetics and Basal Levels of the Transcriptome Response During Effector-Triggered Immunity in Arabidopsis are mainly controlled by Four Immune Signaling Sectors J Bioinform Syst Biol 6 2023 347 363
Hillmer R.A. Tsuda K. Rallapalli G. Asai S. Truman W. Papke M.D. Sakakibara H. Jones J.D.G. Myers C.L. Katagiri F. The highly buffered Arabidopsis immune signaling network conceals the functions of its components PLoS Genet. 13 2017 e1006639 10.1371/journal.pgen.1006639
Huffaker A. Pearce G. Ryan C.A. An endogenous peptide signal in Arabidopsis activates components of the innate immune response Proc. Natl. Acad. Sci. USA 103 2006 10098 10103 10.1073/pnas.0603727103 16785434
Jabs T. Dietrich R.A. Dangl J.L. Initiation of runaway cell death in an Arabidopsis mutant by extracellular superoxide Science 273 1996 1853 1856 10.1126/science.273.5283.1853 8791589
Jones J.D.G. Dangl J.L. The plant immune system Nature 444 2006 323 329 10.1038/nature05286 17108957
Kermack W.O. McKendrick A.G. A contribution to the mathematical theory of epidemics Proceedings of the Royal Society of London. Series A 115 1927 700 721
Khan M. Seto D. Subramaniam R. Desveaux D. Oh, the places they'll go! A survey of phytopathogen effectors and their host targets Plant J. 93 2018 651 663 10.1111/tpj.13780 29160935
Király L. Künstler A. Höller K. Fattinger M. Juhász C. Müller M. Gullner G. Zechmann B. Sulfate supply influences compartment specific glutathione metabolism and confers enhanced resistance to Tobacco mosaic virus during a hypersensitive response Plant Physiol. Biochem. 59 2012 44 54 10.1016/j.plaphy.2011.10.020 22122784
Krol E. Mentzel T. Chinchilla D. Boller T. Felix G. Kemmerling B. Postel S. Arents M. Jeworutzki E. Al-Rasheid K.A.S. Perception of the Arabidopsis danger signal peptide 1 involves the pattern recognition receptor AtPEPR1 and its close homologue AtPEPR2 J. Biol. Chem. 285 2010 13471 13479 10.1074/jbc.M109.097394 20200150
Liu X. Igarashi D. Hillmer R.A. Stoddard T. Tsuda K. Myers C.L. Katagiri F. Dynamic decomposition of transcriptome responses during plant effector-triggered immunity revealed conserved responses in two distinct cell populations Preprint at bioRxiv 2022 10.1101/2022.1112.1130.522333
Liu Z. Wu Y. Yang F. Zhang Y. Chen S. Xie Q. Tian X. Zhou J.M. BIK1 interacts with PEPRs to mediate ethylene-induced immunity Proc. Natl. Acad. Sci. USA 110 2013 6205 6210 10.1073/pnas.1215543110 23431184
Lu Y. Tsuda K. Intimate Association of PRR- and NLR-Mediated Signaling in Plant Immunity Mol. Plant Microbe Interact. 34 2021 3 14 10.1094/MPMI-08-20-0239-IA 33048599
Mackey D. Holt B.F. Wiig A. Dangl J.L. RIN4 interacts with Pseudomonas syringae type III effector molecules and is required for RPM1-mediated resistance in Arabidopsis Cell 108 2002 743 754 10.1016/s0092-8674(02)00661-x 11955429
Mackey D. Belkhadir Y. Alonso J.M. Ecker J.R. Dangl J.L. Arabidopsis RIN4 is a target of the type III virulence effector AvrRpt2 and modulates RPS2-mediated resistance Cell 112 2003 379 389 10.1016/s0092-8674(03)00040-0 12581527
Mathew I.E. Agarwal P. May the Fittest Protein Evolve: Favoring the Plant-Specific Origin and Expansion of NAC Transcription Factors Bioessays 40 2018 e1800018 10.1002/bies.201800018
Mine A. Seyfferth C. Kracher B. Berens M.L. Becker D. Tsuda K. The Defense Phytohormone Signaling Network Enables Rapid, High-Amplitude Transcriptional Reprogramming during Effector-Triggered Immunity Plant Cell 30 2018 1199 1219 10.1105/tpc.17.00970 29794063
Nobori T. Monell A. Lee T.A. Zhou J. Nery J. Ecker J.R. Time-resolved single-cell and spatial gene regulatory atlas of plants under pathogen attack Preprint at bioRxiv 2023 10.1101/2023.1104.1110.536170
O'Malley R.C. Huang S.S.C. Song L. Lewsey M.G. Bartlett A. Nery J.R. Galli M. Gallavotti A. Ecker J.R. Cistrome and Epicistrome Features Shape the Regulatory DNA Landscape Cell 166 2016 1598 10.1016/j.cell.2016.08.063 27610578
Prabakaran R. Jemimah S. Rawat P. Sharma D. Gromiha M.M. A novel hybrid SEIQR model incorporating the effect of quarantine and lockdown regulations for COVID-19 Sci. Rep. 11 2021 24073 10.1038/s41598-021-03436-z
Ramsay J.O. Silverman B.W. Functional Data Analysis 2005 Springer
Rescigno A. Synthesis of a Multicompartmented Biological Model Biochim. Biophys. Acta 37 1960 463 468 14437430
R_Core_Team A Language and Environment for Statistical Computing 2021 R Foundation for Statistical Computing
Sorenson R.S. Deshotel M.J. Johnson K. Adler F.R. Sieburth L.E. mRNA decay landscape arises from specialized RNA decay substrates, decapping-mediated feedback, and redundancy Proc. Natl. Acad. Sci. USA 115 2018 E1485 E1494 10.1073/pnas.1712312115 29386391
Su J. Yang L. Zhu Q. Wu H. He Y. Liu Y. Xu J. Jiang D. Zhang S. Active photosynthetic inhibition mediated by MPK3/MPK6 is critical to effector-triggered immunity PLoS Biol. 16 2018 e2004122 10.1371/journal.pbio.2004122
Tsuda K. Sato M. Stoddard T. Glazebrook J. Katagiri F. Network properties of robust immunity in plants PLoS Genet. 5 2009 e1000772
Tsuda K. Qi Y. Nguyen L.V. Bethke G. Tsuda Y. Glazebrook J. Katagiri F. An efficient Agrobacterium-mediated transient transformation of Arabidopsis Plant J. 69 2012 713 719 10.1111/j.1365-313X.2011.04819.x 22004025
van der Hoorn R.A.L. Kamoun S. From Guard to Decoy: a new model for perception of plant pathogen effectors Plant Cell 20 2008 2009 2017 10.1105/tpc.108.060194 18723576
Xin X.F. Nomura K. Aung K. Velásquez A.C. Yao J. Boutrot F. Chang J.H. Zipfel C. He S.Y. Bacteria establish an aqueous living space in plants crucial for virulence Nature 539 2016 524 529 10.1038/nature20166 27882964
Yamaguchi Y. Huffaker A. Bryan A.C. Tax F.E. Ryan C.A. PEPR2 is a second receptor for the Pep1 and Pep2 peptides and contributes to defense responses in Arabidopsis Plant Cell 22 2010 508 522 10.1105/tpc.109.068874 20179141
Zaghlool A. Niazi A. Björklund Å.K. Westholm J.O. Ameur A. Feuk L. Characterization of the nuclear and cytosolic transcriptomes in human brain tissue reveals new insights into the subcellular distribution of RNA transcripts Sci. Rep. 11 2021 4076 10.1038/s41598-021-83541-1 33603054
Zhu J. Lolle S. Tang A. Guel B. Kvitko B. Cole B. Coaker G. Single-cell profiling of Arabidopsis leaves to Pseudomonas syringae infection Cell Rep. 42 2023 112676 10.1016/j.celrep.2023.112676
