The kinetics and basal levels of the transcriptome response during Effector-Triggered Immunity in Arabidopsis are mainly controlled by four immune signaling sectors

To observe the transcriptome response during Effector-Triggered Immunity (ETI) without complications from any other pathogen factors or heterogeneously responding cell populations, we transgenically and conditionally expressed the Pseudomonas syringae effector AvrRpt2 in Arabidopsis leaves. We studied this ETI-specific, cell-autonomous transcriptome response in 16 exhaustively combinatorial genetic backgrounds for the jasmonate (JA), ethylene (ET), PAD4, and salicylate (SA) immune signaling sectors. Removal of some or all four sectors had relatively small impacts on the intensity of the overall ETI transcriptome response (1972 upregulated and 1290 downregulated genes). Yet, we found that the four signaling sectors strongly affect the kinetics of the ETI transcriptome response based on analysis of individual genes via time-course modeling and of the collective behaviors of the genes via a PCA-based method: the PAD4 sector alone and the JA;SA sector interaction (defined by the averaging model) accelerated the response, while the ET;SA sector interaction delayed it. The response acceleration by the PAD4 sector or the JA;SA sector interaction was consistent with their positive contributions to ETI measured by pathogen growth inhibition. The responsive genes overlapping between ETI and Pattern-Triggered Immunity (PTI) had distinct regulatory trends regarding the four sectors, indicating different regulatory circuits in upstream parts of ETI and PTI signaling. The basal mRNA levels of most ETI-upregulated genes, but not downregulated genes, were predominantly positively regulated by the PAD4;SA sector interaction. This detailed mechanistic decomposition of the roles of four signaling sectors allowed us to propose a potential regulatory network involved in ETI signaling.


SUMMARY
To observe the transcriptome response during Effector-Triggered Immunity (ETI) without complications from any other pathogen factors or heterogeneously responding cell populations, we transgenically and conditionally expressed the Pseudomonas syringae effector AvrRpt2 in Arabidopsis leaves. We studied this ETI-specific, cell-autonomous transcriptome response in 16 exhaustively combinatorial genetic backgrounds for the jasmonate (JA), ethylene (ET), PAD4, and salicylate (SA) immune signaling sectors. Removal of some or all four sectors had relatively small impacts on the intensity of the overall ETI transcriptome response (1972 upregulated and 1290 downregulated genes). Yet, we found that the four signaling sectors strongly affect the kinetics of the ETI transcriptome response based on analysis of individual genes via time-course modeling and of the collective behaviors of the genes via a PCA-based method: the PAD4 sector alone and the JA;SA sector interaction (defined by the averaging model) accelerated the response, while the ET;SA sector interaction delayed it. The response acceleration by the PAD4 sector or the JA;SA sector interaction was consistent with their positive contributions to ETI measured by pathogen growth inhibition. The responsive genes overlapping between ETI and Pattern-Triggered Immunity (PTI) had distinct regulatory trends regarding the four sectors, indicating different regulatory circuits in upstream parts of ETI and PTI signaling. The basal mRNA levels of most ETI-upregulated genes, but not downregulated genes, were predominantly positively regulated by the PAD4;SA sector interaction. This detailed mechanistic decomposition of the roles of four signaling sectors allowed us to propose a potential regulatory network involved in ETI signaling.

INTRODUCTION
Regulatory mechanisms of plant immunity during interactions with pathogens are complex. This is because plant immune systems must have a high level of resilience to withstand assaults from fast-evolving pathogens (Katagiri, 2018). A common approach to a complex problem is reductionism: reducing a complex problem to a combination of a small number of simpler problems. In plant immunity, a common reduced conceptual framework is a sequence of Pattern-Triggered Immunity (PTI), Effector-Triggered Susceptibility (ETS), and Effector-Triggered Immunity (ETI) (Dodds & Rathjen, 2010;Jones & Dangl, 2006). PTI is elicited when pattern recognition receptors (PRRs) on the plant cell membrane recognize the cognate molecular patterns, which are conserved among phylogenetically similar microbes or are plant molecules associated with damage. For example, a part of bacterial flagellin, flg22, is recognized as a microbe-associated molecular pattern (MAMP) by the Arabidopsis PRR FLS2 (Gómez-Gómez et al, 1999;Zipfel et al, 2004). Pathogens well-adapted to the host plant deliver effector molecules into the plant cell and compromise PTI signaling, which is the state of ETS. Plants have intracellular ETI receptors called resistance (R) proteins that sometimes could directly or indirectly recognize the cognate effectors delivered. For example, the Pseudomonas syringae effector AvrRpt2 is recognized by the Arabidopsis ETI receptor RPS2 (Bent et al, 1994;Mindrinos et al, 1994). This recognition leads to elicitation of ETI.
To apply reductionism to ETI according to this framework, the state of ETI needs to be studied in isolation from PTI and ETS. If an ETI-eliciting effector is delivered from a pathogen, the pathogen necessarily presents MAMPs, and consequently PTI is elicited as well. Since ETI and PTI interact in a highly context dependent manner (Hatsugai et al, 2017;Ngou et al, 2021;Wang et al, 2023;Yuan et al, 2021), subtraction of the quantitative measure of plant immunity with a pathogen strain (a mixture of PTI and ETS states) from that with the congenic pathogen strain carrying an ETI-eliciting effector (a mixture of PTI, ETS, and ETI states) does not provide a measure of the ETI-only state. Another complexity for studying plant immunity is heterogenous states among plant cells. During ETI, the plant cells that directly recognize an ETI-eliciting effector typically undergo the hypersensitive response (HR), which is a cell-autonomous programmed cell death (Coll et al, 2011). The surrounding plant cells activate a later immune 4 response based on some secondary signal(s) from the autonomously responding plant cells and/or diffusible pathogen factors (Liu et al, 2022;Lu & Tsuda, 2021). Spatially separating the responses of these two cell populations experimentally via single cell RNA-seq with good time resolution and adequate replication is technically and financially challenging. Artificial expression of an ETI-eliciting effector in plant cells circumvents these issues (Hatsugai et al., 2017;Ngou et al., 2021;Yuan et al., 2021): the ETI-only, cell-autonomous-response-only state can be achieved. In a whole-plant application, use of a chemical-inducible promoter controlling expression of an effector is a practical method, e.g., (McNellis et al, 1998;Tsuda et al, 2012).
Another area in which reductionism can be applied to the study of ETI is dissection of the ETI signaling network. Signaling machineries involved in the network interact with each other, making the behavior of the signaling system complex: particularly, backup mechanisms (buffering interactions) in the network provide resilience against disabling attack on some of the machineries (Hillmer et al, 2017;Katagiri, 2018;Tsuda et al, 2009). We have reduced the signaling network to a network of five major immune signaling sectors, the jasmonate (JA), ethylene (ET), PAD4, and salicylate (SA) sectors and the ETI-Mediating, PTI-Inhibited sector (EMPIS), in Arabidopsis (Hatsugai et al., 2017;Hillmer et al., 2017;Tsuda et al., 2009).
Thus, these three machineries were largely sufficient for ETI signaling in the five-sector context. An important step in reductionism is to synthesize understanding of the original complex problem by combining what we learned from several simpler problems. However, when no generic rules about how to combine them are known, it is necessary to measure what happens when they are combined in every combination, i.e., study their interactions in an exhaustive manner. Regarding four signaling sectors, the JA, ET, PAD4, and SA sectors, we have been using 16 exhaustively combinatorial genotypes among the hub genes of the four sectors to estimate all possible interactions among them (Hillmer et al., 2017;Tsuda et al., 2009). This 5 essentially reconstitutes the intact network sector by sector, so we call this analytical approach network reconstitution (formerly, signaling allocation) (Katagiri, 2017). We improved a general linear model describing network reconstitution so that it allows an intuitive and consistent interpretation of higher order interactions among sectors (Katagiri, 2022). The new model is called the averaging model, and we call the analysis Network Reconstitution via Averaging Model (NRAM). The definition of interactions in the averaging model is different from that in an additive model, so a ";" is used to indicate an averaging model interaction, e.g., JA;SA interaction.
Here, we collected Arabidopsis RNA-seq data in response to Estradiol (Ed) -inducible in planta expression of AvrRpt2 in 16 exhaustively combinatorial genetic backgrounds for the four major immune signaling sectors, the JA, ET, PAD4, and SA sectors. We applied NRAM with the four sectors to the cell-autonomous ETI transcriptome response. The transcriptome response was characterized by the kinetic parameters of the peak amplitude and time of the response timecourse for each of 1972 ETI-upregulated and 1290 ETI-downregulated genes. The four sectors did not strongly affect the overall response intensity, consistent with the fact that EMPIS (Hatsugai et al., 2017), which could be inhibited by PTI, was not disabled. However, the PAD4 sector alone or the JA;SA sector interaction, which positively contribute to AvrRpt2-ETI, accelerated the responses of both upregulated and downregulated genes. This observation suggests the importance of a fast response during AvrRpt2-ETI for strong immunity. In contrast, the ET;SA sector interaction delayed these responses. We applied similar analyses to a previously published flg22-PTI transcriptome response data set (Hillmer et al., 2017). Although the responsive genes were highly overlapping between AvrRpt2-ETI and flg22-PTI, regulation of the transcriptome response via the four signaling sectors was distinct between AvrRpt2-ETI and flg22-PTI. We also found that the basal mRNA levels of ETI-upregulated, but not ETIdownregulated, genes were predominantly and positively regulated by the PAD4;SA sector interaction. We propose a regulatory network model that can explain our observations regarding ETI transcriptome response regulation across the four signaling sectors.

Genotype nomenclatures
For simplicity, we use the following genotype nomenclature. The names of transgenes for Edinduced expression are preceded by Ed-, such as Ed-AvrRpt2. All genotypes carried the Ed-AvrRpt2 transgene, except the GUS genotype carrying Ed-GUS in a wild-type genetic background. The wild-type alleles regarding of the hub genes of the four signaling sectors, the JA, ET, PAD4, and SA sectors, are shown by single letters, J, E, P, and S, respectively. The corresponding mutant alleles are shown by the letter x for the corresponding position. For example, JEPS represents the wild type alleles for all four genes, and JxPx represents the wildtype alleles for the hub genes of the JA and PAD4 sectors and mutant alleles for those of the ET and SA sectors. This genotype nomenclature for the four sectors, using "x", makes it clear which sectors remain functional in a particular genotype, which is important in mechanistic interpretations of the sectors (Katagiri, 2022;Tsuda et al., 2009). The r2r1 genotype carried mutant alleles of the R genes, rps2 and rpm1, while its four sector hub genes had wild-type alleles.

Overview of the RNA-seq data
The RNA-seq data for each genotype * treatment * time combination were consistent across biological replicates by visual inspection, indicating high reproducibility of the data ( Figure S1).
The mean estimates for four genotypes, JEPS, xxxx, r2r1, and GUS across the treatment and the 7 time for 18970 well-expressed genes were subjected to PCA. The results are shown in a PC1 x PC2 plot ( Figure 1A). While the PC1 values of JEPS and xxxx with Ed showed large variation in a time-dependent manner, the PC1 values among the negative controls (r2r1, and GUS with Ed and all mocks) showed only small variation. Therefore, PC1, which captures 45.1% of the data variance, mostly represents the ETI-specific transcriptome response. On the other hand, PC2, which represents 17.3% of the variance, mainly captures circadian/diurnal responses (Yang et al, 2020) in the negative controls (Text S1 for details). Although the ETI transcriptome response in xxxx was very similar to the response in JEPS as the trajectories along their time points were very close in the PC1 x PC2 space, the xxxx response was generally slower than the JEPS response (compare "E3" to "E6" between the genotypes). These overview observations regarding the ETI-response in JEPS and xxxx are in agreement with previous observations when AvrRpt2 was delivered from P. syringae (Mine et al, 2018).

AvrRpt2 mRNA level variation did not cause the ETI transcriptome response variation across the combinatorial genotypes
Since AvrRpt2 is the signal that initiates an ETI response in our experimental system, if AvrRpt2 expression was significantly leaky at 0 hpt or affected by the combinatorial genotypes, it could complicate interpretations of the results. Thus, before further analyzing the ETI transcriptome response across the combinatorial genotypes, we examined these possibilities.
First, uninduced, leaky expression of AvrRpt2 was investigated because leaky AvrRpt2 expression in the Ed-AvrRpt2-carrying lines could elicit a constitutive, weak ETI response.
Although our RNA-seq read depth per library was generally limited (1.2 to 8.4 million aligned reads per library; median 3.4 million), 96 out of 102 libraries for the genotypes carrying Ed-AvrRpt2 at 0 hpt had no read counts for AvrRpt2, indicating very low leaky expression of AvrRpt2 in general. Then, we compared mRNA level differences gene by gene between every pair of JEPS, r2r1, and GUS genotypes at 0 hpt. Figure 1B shows a heatmap for three pairwise comparisons with 158 genes that showed significant differences in at least one of the pairwise comparisons. If leaky AvrRpt2 expression caused a weak ETI response, the genes for a weak 8 ETI response should show little difference in r2r1/GUS and significant and similar differences in r2r1/JEPS and GUS/JEPS since a weak ETI response would be expected only in JEPS. There was no gene that fit this profile of a potential weak ETI response. All genotypes except GUS shared the same Ed-AvrRpt2 transgene locus introduced by genetic crosses from a single transgenic line (see Methods). Thus, we conclude that there was no appreciable effect of leaky AvrRpt2 expression in our lines.
Second, we investigated induced AvrRpt2 mRNA levels across the combinatorial genotypes. The induced AvrRpt2 mRNA levels (the log2(Ed/mock) mean estimates) indeed varied substantially across the genotypes ( Figure 2A). However, the induced AvrRpt2 mRNA level was generally strongly anticorrelated with those of ETI-responsive genes across the genotypes ( Figure 2B), indicating that variation in the AvrRpt2 mRNA level was not a major cause of the ETI transcriptome response variation we observed across the genotypes. Figure 3A shows a heatmap of the log2(Ed/mock) mean estimates for 16 combinatorial genotypes at 7 time points (0, 1, 2, 3, 4, 5, and 6 hpt). It contains 1972 ETI-upregulated and 1290 downregulated genes, which were well-modeled for the time course and ordered according to the model-fit peak times in JEPS (see below for the time-course model). The combinatorial genotypes from JEPS to xxxx did not show large differences in the log2(Ed/mock) time courses across the genes. This was not surprising because EMPIS can mediate relatively intact ETI when it is not inhibited by PTI signaling (i.e., the four signaling sectors are not absolutely required for an ETI response when EMPIS is functional) (Hatsugai et al., 2017). However, when the log2(Ed/mock) value for xxxx were subtracted from those for JEPS at each time point, the peaks of the differences were generally earlier than the peaks of the values in xxxx ( Figure 3B), indicating that the peak times were earlier in JEPS than xxxx for most of the upregulated and downregulated genes. Figure 3C shows the PC1 x PC2 plot of genotype * time data resulting from PCA of the values used in Figure 3A. PC1 captured most of the ETI response across the genotypes (83.6% variance). In the PC1 x PC2 plot (together, 90% variance), the response trajectories of the time points for the genotypes were all very similar, monotonously increasing along PC1, and differed in how fast the PC1 value increased. This observation confirms that the ETI-transcriptome response was similar across the genotypes and that the major difference across the genotypes was response kinetics. Subsequently, the PC1 value in Figure 3C is used as the representative value for the ETI-response of each genotype (designated the "ETI-PC1" value). Figure 3D shows the time courses of ETI-PC1. Roughly speaking, the more sectors that were functional, the faster the  Table S1 shows the fitted parameter values and the fitted log2(Ed/mock) values at each time point). The model parameters were the peak amplitude (positive and negative for upregulated and downregulated), the peak time, the peak shape, and the time lag. The log2(Ed/mock) time courses were well-modeled for 1972 upregulated and 1290 downregulated genes. However, since the peak characteristics of the models with peak times well after the last data time point (6 hpt) are not reliable, 866 upregulated and 234 downregulated genes with estimated peak times < 7.5 hpt in at least 13 out of 16 genotypes were selected for NRAM analysis of the peak amplitude and time ( Figures 4B and 4C). The left-most lane in the NRAM heatmap for the peak amplitude shows the log2-ratio of the peak amplitudes between JEPS and xxxx for each gene ("PA,JEPS/xxxx" in Figure 4B). The general trend was that the absolute value of the peak amplitude was larger in JEPS than xxxx: log2(JEPS/xxxx) medians 0.51 and -0.74 for the upregulated and downregulated genes, respectively. In the case of the peak time differences between JEPS and xxxx ("PT,JEPS-xxxx" in Figure 4C), in agreement with the observations in "JEPS/xxxx" in Figure 3B (the same lanes are also shown in Figure 4B as a reference), the general trend was delays in the peak time in xxxx compared to JEPS: medians 0.92-and 0.77-hour delays for the upregulated and downregulated genes, respectively. These values, "PA,JEPS/xxxx" and "PT,JEPS-xxxx", were subjected to NRAM-decomposition for the peak amplitude and time, respectively. Since the number of downregulated genes analyzed was These trends observed in NRAM analysis were only for the well-modeled upregulated, earlypeak genes, which are about half of the well-modeled upregulated genes, and no trends were evident for the small number of well-modeled downregulated, early-peak genes. Since the ETI-PC1 captured most of the ETI response and its timing ( Figure 3C), we subjected the ETI-PC1 values at 4 hpt to NRAM analysis. The time of 4 hpt was chosen because most genes started responding by this time, and most genes had not reached their peaks by this time, so that this time point likely captures the response kinetics of the gene sets with little bias. In fact, the obtained ETI-PC1 NRAM profiles at 4 hpt for 1972 ETI-upregulated genes ( Figure 4F) and 1290 downregulated genes ( Figure 4G) were very similar to the peak time NRAM mean profiles of 866 upregulated genes ( Figure 4E). The similar mean trends of the NRAM profiles between the upregulated and downregulated genes indicate that the peak times of the upregulated and downregulated genes were on average regulated essentially in the same manner by the four signaling sectors. Using bootstrapping, the 95% confidence intervals of the NRAM terms were estimated for the upregulated and downregulated genes. The small 95% confidence intervals relative to the mean values suggest that the observed NRAM profile trends are very common across the genes in each of the upregulated and downregulated gene sets, i.e., the peak times of most of the upregulated and downregulated genes were regulated essentially in the same manner by the four signaling sectors.

The four signaling sectors mainly affect the response kinetics
The roles of the four sectors in regulation of PTI-responsive genes were clearly different from those in regulation of ETI-responsive genes.
A similar analysis was conducted with the flg22-PTI transcriptome response (data from (Hillmer et al., 2017)) for the genes in common between this response and the 3262 ETI-responsive wellmodeled genes (602 ETI-upregulated and 594 downregulated genes were common; Figure 5).
The gamma-pdf time course model (with no time lag parameter) was fit to the estimated log2(flg22/mock) data at 2, 3, 5, and 9 hpt for each of the 16 combinatorial genotypes for each gene ( Figures S4 and S5 show the model fit by heatmaps and time courses, respectively; Table   S2 shows the fitted parameter values and the fitted log2(flg22/mock) values at each time point).
The data at 1 hpt were not used because the response at 1 hpt likely represents a general stress response rather than an immune-specific response (Bjornson et al, 2021). The data at 18 hpt were not used because we wanted to compare the responses between ETI and PTI over similar time ranges. The genes upregulated or downregulated in ETI and PTI were generally consistent across the common genes ("xxxx" in Figure 5A). The collective effect of the four sectors on the peak amplitude was limited for most genes (faint colors in "PA,JEPS/xxxx" compared to "xxxx" in Figure 5A), and the collective effect was not consistent across the upregulated or downregulated genes (reddish and blueish color bands are mixed in "PA,JEPS/xxx"). Furthermore, no consistent pattern was evident in the NRAM results for the PTI peak amplitude. On the other hand, consistent patterns were evident in the NRAM results for the PTI peak time: The ET sector accelerated and the JA;SA sector interaction delayed the PTI transcriptome response ( Figure   5B), which was confirmed by the mean values of the NRAM terms ( Figure 5C). When PCA was applied to the PTI transcriptome response data, PC1 appeared not to capture the timing of the response (the genotype-specific trajectories were not highly overlapping; Figure 5D). Therefore, the PC1 values cannot be used for analysis of the peak time in the PTI response. In short, control of the peak amplitude and time by the four signaling sectors is distinct between the ETI and PTI transcriptome responses.

12
The basal mRNA levels of most ETI-upregulated genes were largely positively controlled by the PAD4;SA sector interaction We also investigated the influence of the four signaling sectors on the uninduced, basal mRNA levels of the ETI-responsive genes ( Figure 6). The basal mRNA levels of most ETI-upregulated genes but not the downregulated genes were positively regulated by the PAD4;SA sector interaction ( Figure 6A). Thus, the basal levels of the ETI-upregulated genes are likely positively regulated by the PAD4;SA sector interactions. We examined the converse of this observation, whether the genes whose basal mRNA levels are positively regulated by the PAD4;SA sector interactions are likely to be ETI-upregulated genes. When the 4562 genes whose basal mRNA levels were affected in some way by the four signaling sectors were subjected to NRAM, the genes with positive PAD4;SA sector interaction values were strongly associated with the ETIupregulated genes ( Figure 6B). Thus, the gene set for the ETI-upregulated genes and the set of genes whose basal mRNA levels are positively regulated by the PAD4;SA sector interactions are largely overlapping.

DISCUSSION
In this study, we combined simplification of an ETI experimental system by in planta expression of the ETI-eliciting effector AvrRpt2, dynamical analysis of the ETI transcriptome response, and decomposition of the signaling sector contributions to the transcriptome response by NRAM.
This multifaceted analytical design allowed us to reveal how the dynamical characteristics of the ETI-transcriptome responses are regulated through specific signaling sectors and their interactions. Similar analysis of the PTI transcriptome response indicated that dynamical characteristics of the overlapping gene set were regulated differently from those in ETI. Furthermore, the ETI upregulated genes were characterized by positive regulation of their basal mRNA levels via the PAD4;SA sector interaction.
The wild type alleles of RPS2 and/or RPM1 altered the transcriptome without ETI elicitation ( Figure 1B). The genes downregulated by RPS2 and/or RPM1 (R2R1_down) overlap with ETIupregulated genes (ETI_up; Figure 7A). This tendency is also detectable in Figure 1A, which shows that the data points for the r2r1 genotype have slightly higher PC1 values than the other 13 data points for no ETI elicitation in the genotypes containing the R2R1 wild-type alleles. Since PC1 mainly represents the ETI response, the slightly higher PC1 values in r2r1 correlate with downregulation of the ETI response by RPS2 and/or RPM1. Since it has been reported that functional alleles of RPS2 downregulate the mRNA levels of some defense-related genes (MacQueen et al, 2016), this observation could be explained solely by the RPS2 effect. We only saw 12 genes significantly upregulated by RPS2 and/or RPM1, and they have virtually no overlap with ETI-responsive genes ( Figures 7A and 7B). The R2R1 downregulated gene set (142 genes) was mostly included in a larger set of genes whose basal mRNA levels were positively regulated by the PAD4;SA interaction (738 genes; Figure 7C), suggesting that RPS2 and/or RPM1 lowers the basal mRNA levels of the R2R1 downregulated genes by weakly interfering with the interaction between the PAD4 and SA sectors.
The four-sector influence on AvrRpt2 mRNA level upregulation is completely different from those in most ETI-upregulated genes ( Figure 2B). The genotype difference in the AvrRpt2 mRNA level was already clear at 1 hpt ( Figure 2A) while not much ETI transcriptome response was observed before 2 hpt (e.g., Figures 1A, 3A, and 3C). These two observations suggest that the genotype difference in the AvrRpt2 mRNA level is not a consequence of a difference in ETI response across the genotypes. The difference in the XVE protein level (the artificial TF for the Ed-inducible promoter activation; (Zuo et al, 2000)) likely explains it as the mean XVE mRNA level ( Figure 2C) positively correlates with the AvrRpt2 mRNA accumulation rate ( Figure 2D). The ETI signaling system is built in a way that allows three signaling mechanisms, the EMPIS only, the PAD4 sector only, and the JA and SA sectors together, to back up each othernone of the three is essential, but any one of the three is sufficient (Hatsugai et al., 2017;Katagiri, 2022;Tsuda et al., 2009). In our experimental system, EMPIS was not disabled, and thus the influence of the four signaling sectors we observed was limited to modification of the EMPIS-mediated ETI response. Currently the only known way to disable EMPIS is to elicit PTI signaling (Hatsugai et al., 2017). However, this is not compatible with our objective of studying ETI in isolation from PTI. In addition, PTI signaling interferes with inducible promoter systems (Igarashi et al, 2013), so it would be experimentally difficult to exclude PTI influence on Ed-AvrRpt2 expression. If we knew how to disable EMPIS independent of PTI, such as by a genetic 14 mutation, it would be interesting to investigate how the ETI response is restored from a completely disabled state through each of the three ways.
The peak amplitude and time in the ETI transcriptome response were regulated differently by the four signaling sectors. We estimated the four signaling sector influence on the peak amplitude and time for 866 upregulated genes with relatively early peak times ( Figures 4B-E). For most of these genes, the peak amplitude was regulated by the PAD4 and SA sectors positively and by the PAD4;SA sector interaction negatively. For the same genes, the peak time was accelerated by the PAD4 sector and the JA;SA sector interaction. Using the ETI-PC1 value, we extended the observation of the four signaling sector influence on the peak time to most of the upregulated and downregulated genes. The PAD4 sector and the JA;SA sector interaction significantly and positively contribute to ETI-mediated resistance against P. syringae (Katagiri, 2022;Tsuda et al., 2009): The influences of the two mechanisms were obvious in the case of AvrRpt2-ETI, in which EMPIS was inhibited by PTI signaling. They were also evident even when EMPIS was active in the case of AvrRpm1-ETI ( Figures 7A and 7B in (Katagiri, 2022). The latter suggests that acceleration of the ETI response is important to enhance ETI-mediated resistance against P.

syringae.
One simple signaling circuit that allows regulation of the peak amplitude and time separately is a feedforward loop with a linear combination at the signal convergence point ( Figure 8A). Using a multi-compartment model (Godfrey, 1983;Kermack & McKendrick, 1927;Liu et al., 2022;Prabakaran et al, 2021;Rescigno, 1960) for a feedforward loop, the peak time can be changed by changing the linear coefficients and the decay rate of the converging node ( Figure 8B). The peak amplitude can be changed by these three model parameters and the input amplification ratio parameter before the loop. Four model parameters can be used to determine two model characteristics of the peak amplitude and time to change, and thus the feedforward model is highly flexible in determining the peak amplitude and time.
We observed a highly common characteristic that the basal mRNA levels of ETI-upregulated genes were positively regulated by the PAD4;SA sector interactions ( Figure 6). This characteristic was not shared by ETI-downregulated genes, strongly suggesting that this basal 15 level regulation acts after signals for upregulation and downregulation are split. Such a basal level regulation could be explained by an additional non-ETI-responding input or upregulation of the input amplification ratio parameter shared among the upregulated genes ( Figure 8C, top and bottom, respectively). The additional input does not change the difference between the peak level and the basal level on a linear scale ( Figure 8D, top left). The input amplification ratio can multiply the basal level and the peak level equally. In the latter scenario, the difference between the induced and basal levels is unchanged on a log scale ( Figure 8D, bottom right), which indicates that the fold change (induced / basal) is not changed when the input amplification ratio is changed. The data in this study does not have sufficient statistical power to test these two possible scenarios.
In defense priming situations in Arabidopsis, the mRNA levels of defense-related genes after mock or pathogen inoculation in secondary leaf tissue were compared between prior mock and a defense primer (Pto DC3000 or pipecolic acid) treatments of primary leaf tissue (Bernsdorff et al, 2016). Defense primer treatments induce the systemic acquired resistance (SAR) state in the secondary tissue. Generally, SAR induction changed the basal level but did not strongly change the fold change of pathogen-upregulated genes in the secondary tissue, which is consistent with basal level regulation by the input amplification ratio (Figures 6 and 8 in (Bernsdorff et al., 2016)): defense priming substantially increases the input amplification ratio. Since the genes upregulated in different types of immunity overlap, it is likely that the basal level regulation we observed with the ETI-upregulated genes is also mediated by regulation of the input amplification ratio.
The transcriptome response regulation by the four signaling sectors was very different between ETI and PTI while the overlap of upregulated gene sets was large. We proposed a WRKY network activation hypothesis (Liu et al., 2022). In this hypothesis, PTI and ETI signaling have separate upstream parts so that they can be regulated separately, and then the signals are fed into different entry points of the transcriptionally interconnected WRKY network. These inputs through different entry points would activate highly overlapping WRKY subgroups although the specific WRKY members that are activated are not the same, explaining the overlapping but not identical upregulated gene sets.
We combined two feedforward loops for ETI and PTI for both upregulated and downregulated genes, basal level regulation by the input amplification ratio for upregulated genes, and WRKY network activation for upregulated genes, to propose one possible regulatory network model for ETI and PTI transcriptome responses (Figure 9). Although it was not clear how extendable was the peak amplitude regulation we observed with ETI-upregulated genes with relatively early peak times, we tentatively assumed that this regulation is common among both upregulated and downregulated genes in this model. One ambiguity we left in the model is whether the input amplification ratio-based basal level regulation is part of the WRKY network or upstream of the WRKY network. Since many genes whose basal mRNA levels are positively regulated by the PAD4;SA sector interaction are upregulated in both ETI and PTI, it is likely part of the WRKY network. If this is the case, probably the input amplification ratios of multiple WRKYs are positively regulated by the PAD4;SA sector interactions.
We emphasize that Figure 9 represents just one simple model that can explain all the observations discussed and that there are other possible models. However, we also emphasize that this is a highly testable mechanistic model as it is characterized by a well-defined network architecture and mechanistic model parameters. This level of mechanistic modeling was enabled by simplification of an ETI system by in planta effector expression, dynamical analysis that characterized the response with highly interpretable, quantitative dynamical parameters, and NRAM analysis that allowed proper estimation of interaction interpretations. This approach, comprised of a simplified experimental system, dynamical analysis, and NRAM on dynamical parameters, forms a platform highly applicable to mechanistic studies of other transient response systems that involve interactions of multiple components for their regulation.

Plant materials and growth conditions
An Arabidopsis Col-0 line carrying Ed-AvrRpt2 as a single copy insertion (JEPS; (Tsuda et al., 2012)) was crossed to a dde2-2 ein2-1 pad4-1 sid2-2 quadruple mutant line with Col-0 background (Tsuda et al., 2012). DDE2, EIN2, PAD4, and SID2 are the hub genes of the four immune signaling sectors, the JA, ET, PAD4, and SA sectors, respectively (Alonso et al, 1999;Jirage et al, 1999;von Malek et al, 2002;Wildermuth et al, 2001). Approximately 3500 F2 plants from the cross were progressively screened by PCR-based genotyping (Tsuda et al., 2012;Tsuda et al., 2009) for homozygosity of Ed-AvrRpt2 and either wild-type-or mutant-allele homozygosity for SID2, PAD4, DDE2, and EIN2. For three of the desired genotypes, we only obtained plants with one gene that was heterozygous, so F3 plants from such F2 plants were screened for homozygosity. In this way, all 16 combinatorial genotype lines with the Ed-AvrRpt2 transgene were prepared (JEPS to xxxx). GUS (Tsuda et al., 2012) and r1r2 (Hatsugai et al., 2017) plant lines with Col-0 background were also used. Arabidopsis plants were grown in a controlled environment at 22°C with a 12-hour/12-hour light/dark cycle and 75% relative humidity.

Experimental design of plant treatments
A single 18-genotype set was grown in two 15 cm * 15 cm pots (9 plants per pot), and four 18genotype sets were grown in a single flat of 2 * 4 = 8 pots. In a single biological replicate, four and two flats of plants (16 and 8 18-genotype sets) were grown for Ed and mock treatments, respectively. The flats for Ed and mock treatments were separated to avoid cross-contamination of the treatments. For each of the genotype * treatment * time combinations (180 biological samples), two plants from two different flats were pooled. Thus, in a single biological replicate, 14 and 6 18-genotype sets for Ed and mock treatments (7 and 3 time points), respectively, were used. We used the remaining plants (2 18-genoteype sets for each treatment) as backups in case some plants did not grow well. Plants for three biological replicates were grown one week apart.
In addition, the relative growing positions of 18 genotypes across two pots were randomized across the biological replicates.

Plant treatment and collection
Twenty-four-day old plants were sprayed with 50 µM Ed or mock (both contained 0.25% ethanol in water) three hours after dawn (the time the lights switched on in the growth chamber). The 18 aerial parts of plants were harvested at the indicated time (1, 2, 3, 4, 5, or 6 hpt). Unsprayed plants were used for 0 hpt. The aerial parts of two plants were pooled for each genotype * treatment * time * replicate sample. The harvested plant tissue was immediately flash frozen in liquid N2.

RNA-tag-seq and preprocessing of the RNA-tag-seq data
Each of 540 tissue samples was pulverized and RNA was extracted using TRIzol (Thermo Fisher) and RNeasy (Qiagen) according to the manufacturers' instructions. The RNA-tag-seq libraries were generated according to a 3'-tag sequencing method (Rallapalli et al, 2014). Thirtythree 16-plexed and one 12-plexed libraries were sequenced by the University of Minnesota Genomics Center using an Illumina HiSeq 2000 sequencer. The raw sequence data were preprocessed to obtain read-counts-per-gene data for each library according to (Hillmer et al., 2017). The preprocessed data set included 33603 genes that had at least one read count in at least one library. The raw sequence data and the preprocessed read-counts-per-gene data have been deposited in NCBI GEO (accession: GSE196892).

Estimation of the mean and the continuous mRNA level data
From the ETI-response preprocessed data, genes for which the read count was 0 in more than half of the libraries or for which the 15 th highest read count value across the libraries was fewer than 25 were removed, resulting in 18978 well-expressed genes. Pseudocounts were added to each library proportional to the 90 th percentile of the library (One pseudocount was added to the library with 55 read count as the 90 th percentile value) to have at least one read count in every gene in each library. GLM-NB (negative binomial generalized linear model) with the fixed effect with 18 genotypes * (7 Ed-treated times + 3 mock-treated times) = 180 levels was fit to the pseudocount-added data for each gene with the log-ratio of the 90 th -percentile read count of the libraries over 500 as the offset, to obtain 180 mean estimates, their associated standard errors, and their deviance residuals on a log2 scale. Subsequently all the mRNA level values are on a log2 scale. The deviance residuals were added to their mean estimates to generate approximated continuous log2-transformed mRNA level data (the "continuous mRNA level" data). The deviance residuals were used for this purpose because they normally distribute in a log-scale.
The inverse of the square of the standard errors were used as the weights in linear models when they were fit to the continuous mRNA level data since the continuous data do not preserve the count data errors of the sequence read data.
Using the mean estimates and standard errors from GLM-NB, dynamically responding genes were selected that had significantly different mRNA level values (q < 0.05 (by Storey FDR (Storey et al, 2005)) and |mean difference| > 1 (2-fold)) from 0 hpt at two or more consecutive time points later than 1 hpt in at least one of 18 genotypes, resulting in 10833 dynamical genes.
The experimentally manipulated genes, four genes for combinatorial genotypes, DDE2, EIN2, PAD4, and SID2, two R genes, RPS2 and RPM1, and two transgenes, Ed-AvrRpt2 and XEV, were removed from the figures that show the effects of the Ed-treatment and 18 genotypes unless specifically stated. When 16 combinatorial genotypes were compared, the two R genes, RPS2 and RPM1 were kept in the gene set as they were not manipulated within the 16 genotypes.

Modeling mock time courses using the time course of GUS with Ed for each gene and the mean estimates of log2(Ed/mock) at every time point in every genotype
For each of 10833 dynamical genes in each combinatorial genotype, the mock treatment gives no-ETI response control values. However, we only had three time points for the mock treatment.
We modeled mock mRNA levels for each gene at all seven time points using the time course model for GUS treated with Ed. In Figure 1B, the most prominent pattern was little difference in GUS/JEPS and similar and large differences in r2r1/JEPS and r2r1/GUS, indicating substantial effects of the RPS2 and/or RPM1 genes regardless of the transgene. Therefore, we used the data from GUS instead of r2r1 to model the mock mRNA level values at 7 time points. First, a 4 thorder time polynomial regression model was fit to the GUS continuous mRNA level data.
Second the best polynomial model was selected using the step function and was used as the GUS time course model for the gene. Third, a linear model, ~genotype/time, was fit to the mock 20 continuous mRNA level data for 18 genotypes * 3 time points. Fourth, the best model among ~genotype/time, ~genotype, and ~1 was selected using the step function, and the mean value for 0 hpt of mock was estimated using the best model. Fifth, for each genotype, the estimated 0 hpt of mock was used for the intercept of the time course model of GUS with Ed to obtain the "mock time-course model" (Figure S6). The mock time-course model was used to obtain the mean estimates of the mock values and their standard errors for all genotypes at seven time points.
Then, the mock mean estimates were subtracted from the Ed mean estimates from the GLM-NB model to calculate the mean estimates for log2(Ed/mock) for each gene at every time point in every genotype. The associated standard errors were also calculated. ETI-responsive genes were selected using the criteria that the log2(Ed/mock) values for q < 0.05 and |mean| > log23 (i.e., 3fold) at two or more consecutive time points after 2 hpt in at least one genotype, resulting in 3499 genes.

Modeling the dynamics of the ETI-responsive genes
From 3499 ETI-responsive genes, three non-overlapping gene sets were selected, 404 early.peak genes, 972 early.peak2 genes, and 1912 late.peak genes. For the early.peak genes, genes were selected that had the maximum |log2(Ed/mock)| value among 7 time points later than 2 hpt and earlier than 5 hpt in at least one genotype. For the early.peak2 genes, genes were selected that had the maximum |log2(Ed/mock)| value among 7 time points later than 2 hpt and earlier than 6 hpt in at least one genotype, and then the genes overlapping with early.peak genes were removed. The late.peak genes had the maximum |log2(Ed/mock)| value among 7 time points later than 2 hpt in at least one genotype, and then the genes overlapping with early.peak or early.peak2 genes were removed. Based on visual inspections of the time courses, AT3G03190, AT3G14870, AT3G62950, and ATCG00190 were removed from early.peak genes, and AT1G21160, AT2G14610, AT4G23130, AT5G02320, AT5G10140, AT5G41761, and ATCG01120 were removed from early.peak2 genes, as the time course shapes of some genotypes of these genes did not seem fit the gamma-pdf shape (see the next paragraph). The DDE2, PAD4, and SID2 experimentally manipulated genes and 13 plastid-encoded genes were removed from these gene sets. In the end, we had 399 early.peak, 963 early.peak2, and 1900 21 late.peak genes (total 3262 "modeled" genes; 1972 upregulated and 1290 downregulated genes) selected for the time-course model fitting described below.
We assumed that the shape of the log2(Ed/mock) time course was approximated by the shape of gamma-pdf. Thus, the following function ( ) was fit to each of the early.peak, early.peak2, and late.peak genes.
where is the time in hpt, is the peak amplitude on a log2 scale ( > 0 and < 0 for upregulated and downregulated), is the shape parameter ( > 1; the larger the , the sharper the peak shape), is the peak time in hpt; 0 is the time lag in hpt ( ( < 0 ) = 0); is the baseline to make (0)~0 and (∞)~0. This model is called a gamma-pdf time course model. Example time course models are shown in Figure 4A. The model was fit to the approximated continuous log2(Ed/mock) mRNA ratio data using least square. Among the model parameters, , (log2scaled for a better distribution), and were fit for a time course of each genotype for each gene, while 0 and were made common across the genotypes for each gene. As the statistical power of the data for constraining the parameter values was often limited, the parameter ranges in model fitting were ad hoc constrained. This statistical power issue was particularly severe with the genes and the genotypes with peak times later than 6 hpt, which was the last time point in the data. This was the reason that the genes were first divided into three gene sets according to the peak time, and then the model with different parameter ranges was fit separately for different gene sets. See the R script in Dataset S1 for details.

Network reconstitution via averaging model (NRAM) analysis
For NRAM of the peak amplitudes and times of each gene for the four signaling sectors, the genes with peak times before 7.5 hpt in more than 13 out of 16 combinatorial genotypes were selected from the early.peak and early.peak2 genes (866 "early-modeled" genes). This was because late peak time estimates are not reliable. NRAM analysis was conducted according to 22 (Katagiri, 2022) for the peak amplitudes, the peak times, the basal levels of genes, and the ETI-PC1 values. Briefly, NRAM is a general linear model approach to decompose the phenotype values across the 16 exhaustively combinatorial genotypes. They are decomposed into four signaling sector effects and their interactions. In the averaging model, an interaction (which is denoted using a semicolon, ";") was defined as the difference between the phenotype of the genotype with all sectors in the interaction functional and the average of phenotypes of all genotypes with one of the sectors dysfunctional. This interaction definition makes the averaging model interaction mathematically stable, consistent, and interpretable for any number of sectors involved in the interaction.

Reanalysis of the RNA-seq data for PTI response
The previously published RNA-seq data set for the flg22-PTI response (NCBI GEO accession, GSE78735; (Hillmer et al., 2017)) was reanalyzed similarly to the RNA-seq data set for the ETI response in this study. Details of the analysis are described in Text S2. Briefly, 17423 wellexpressed genes, 12537 dynamical genes, and 2946 PTI-responsive genes were progressively selected, and 2484 genes well-modeled by the gamma-pdf time-course model were obtained. See the R script in Dataset S1 for further details.    values. In addition to the shape parameter, the model is parameterized by the peak amplitude (purple arrow), the time lag (brow arrowhead), and the peak time (blue arrowhead) parameters.
The values for the latter three parameters used for these example time course models are 2.5, 1.5 28 hpt, and 5 hpt, respectively. The value for the baseline parameter, B, is considered subtracted, and therefore its value is zero and not shown. B. NRAM analysis regarding the four sectors was applied to the peak amplitudes estimated by gamma-pdf time course models for each gene during ETI. The left column, indicated as "PA,JEPS/xxxx", shows the ratio of the peak amplitudes between JEPS and xxxx. The set and order of the genes are the same as those in Figure 3A   are also based on the PTI data ("PTI-FC"). The gene set and order are the same as those in Figure 3A (1972 ETI-upregulated and1290 ETI-downregulated genes). The log2 color scale bar shown on the right is the same as that in Figure 3B.     Feedforward loop 1, the current data cannot exclude the possibilities that this peak amplitude regulation occurs at later steps (thus, "?" is attached). The output of Feedforward loop 1 is processed separately for upregulated genes and downregulated genes. The amplification ratio ax for the upregulated gene regulation is positively regulated by the PAD4;SA sector interaction, which was observed as upregulation of the basal mRNA level of ETI-upregulated genes. Then the signal is fed into a transcriptionally interconnected WRKY network (Liu et al., 2022), which upregulates the ETI-upregulated genes. The larger pink dashed ellipse around "WRKY network" indicates the possibility that the PAD4;SA sector interaction control of the amplification ratio ax is included within the function of the WRKY network. A signaling event for the PTI response (right side) is initiated by recognition of flg22 by FLS2 (top right). The generated signal is processed through a separate feedforward loop, Feedforward loop 2, which can regulate the signal kinetics through modulation of coefficients a3 and a4. The signal for regulation of PTIupregulated genes is fed into the WRKY network through different entry-point WRKYs. With different entry points, the WRKY network can upregulate similar but not identical gene sets (Liu et al., 2022).

SUPPLEMENTAL INFORMATION
Text S1. Detailed interpretations of the PCA results shown in Figure 1A.
Text S2. Descriptions of reanalysis of the PTI transcriptome response data. Figure S1. Heatmap of the mRNA levels of all 18978 well-expressed genes x 540 libraries.
Major column divisions by gray vertical lines (36 major columns) divide each genotype * treatment. The major columns 1 to 18 are mock treated and the major columns 19 to 36 are Ed treated. For each treatment, the genotype order is JEPS,xEPS,JxPS,JExS,JEPx,xxPS,xExS,xEPx,JxxS,JxPx,JExx,xxxS,xxPx,xExx,Jxxx,xxxx,r2r1, and GUS. Within a mock major column, the major column is first divided into 3 time points (0, 2, and 5 hpt) and then the time points are divided into three biological replicates. Within an Ed major column, the major column is first divided into 7 time points (0, 1, 2, 3, 4, 5, and 6 hpt) and then the time points are divided into three biological replicates. Note that since the mock major columns have fewer time points, the width of a minor column (each replicate) is wider than that for the Ed major columns. In this way, the major column widths are the same, and the approximately corresponding times across the mock and Ed major columns can be intuitively compared. The color scale for the log2transformed value is shown at the right. Rainbow colors are used so that relatively small differences in the heatmap are readily visible. Figure S2. The fit of the ETI-response gamma-pdf time-course models by heatmaps. In all three heatmaps, the gene set (1972 upregulated and 1290 downregulated genes) and order for the rows and the genotype/time columns are the same as those in Figure 3A. The color code for the log2 ratio is the same as that in Figure 3B. A. A heatmap of the data for the log2 ratio of Ed/mock, which is the same as Figure 3A.   Figures 3A and S2A and Figure S2B. The gene order is the same as that in these figures, top to bottom. The genotype color code and line type code are the same as those used in Figure 3D. Note that the 1-hpt data points were not used in model fitting although 33 some genes show a transient response at 1 hpt. This is because such very early transient responses likely represent general stress responses and are unlikely to represent immune-specific responses (Bjornson et al., 2021). This figure is provided as a separate file since it is very large. Figure S4. The fit of the PTI-response gamma-pdf time-course models for 2484 well-modeled genes by heatmaps. The heatmaps are divided into two for each panel: top, "Common with ETI", the heatmap of the genes overlapping with ETI-response models (602 ETI-upregulated and 594 downregulated genes; Figure 5A) in the same order but with non-overlapping genes removed; bottom, "PTI only", the rest of the genes for the PTI response models ( data and the model. The color code for the log2 ratio is the same as that in Figure 3B. Note that the 1-hpt data points were not used in model fitting although many genes show a transient response at 1 hpt. This is because such very early transient responses likely represent general stress responses and are unlikely to represent immune-specific responses (Bjornson et al., 2021).
The 18-hpt data points were not used in model fitting either, as we tried to compare PTI and ETI responses over a comparable time range. However, a substantial number of genes showed late second-peak responses while the time-course model used can only model single-peak responses.
This is the reason for non-negligible discrepancies at 18 hpt in C. Figure S5. The fit of the PTI-response gamma-pdf time-course models by a time course plot for each of 2484 well-modeled genes. For each gene the data and the model plots are shown side by side (plots for two genes per line). The data and the model values are the same as those used in Figure S4A and Figure S4B. The time axis (the x-axis) of each plot is squire-root-scaled to have the peak shape close to symmetric. The gene order is the same as that in Figure S4A, top to 34 bottom. The genotype color code and the line type code are the same as those used in Figure 3D.
The 1-and 18-hpt data were not used for modeling as explained in the Figure S4 legend. This figure is provided as a separate file since it is very large.  Table S1. The fitted parameter values of the ETI-response gamma-pdf time-course models and the fitted values of the model at 0, 1, 2, 3, 4, 5, and 6 hpt for each genotype for each of 3262 well-modeled genes. The gene order in the row is the same as the gene order in Figure 3A.
Generally, the column names are the genotype name followed by the parameter name or the time point, separated by "_", except the parameters 0 and , which have single values across the genotypes for each gene. The parameters are indicated by "A", "log2s", "pt", "t0", and "B" for , log 2 , , 0 , and in the model, respectively. For example, "JxPS_log2s" indicates the log 2 value for the genotype JxPS, and "xExx_2h" indicates the 2 hpt modeled value for the genotype xExx.