This is the accepted version of the journal article: Salguero Linares, Jose; Serrano, Irene; Ruiz Solaní, Nerea; [et al.]. «Robust transcriptional indicators of immune cell death revealed by spatiotemporal transcriptome analyses». Molecular Plant, Vol. 15, Issue 6 (June 2022), p. 1059-1075. DOI 10.1016/j.molp.2022.04.010 This version is available at https://ddd.uab.cat/record/259326 under the terms of the license Manuscript Click here to access/download;Manuscript;Salguero et al. FINAL.docx Click here to view linked References Robust transcriptional indicators of immune cell death revealed by spatio- temporal transcriptome analyses Running title: Plant HR indicators Short summary: Our work provides a detailed overview of the spatio-temporal transcriptomic landscape of the hypersensitive response (HR), a form of plant-specific immune cell death. Further, a set of robust transcriptional marker genes of the cells that will undergo HR is provided. Jose Salguero-Linaresa,#, Irene Serranob, ,#,†, Nerea Ruiz-Solania, Marta Salas-Gómeza, Ujjal Jyoti Phukana, Victor Manuel Gonzáleza, Martí Bernardo-Fauraa, Marc Vallsa,b, David Rengelb,c,¥,§,*, Nuria S. Colla,d,§,* a Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Campus UAB, Bellaterra, Barcelona, 08193, Spain b LIPM, Universite de Toulouse, INRA, CNRS, 84195 Castanet-Tolosan, France c INRAE, GeT-PlaGe, Genotoul, 31326 Castanet-Tolosan, France (doi: 10.15454/1.5572370921303193E12) d Department of Genetics, Universitat de Barcelona, 08028 Barcelona, Spain e Consejo Superior de Investigaciones Científicas (CSIC), Barcelona, Spain #,§ Equal contributions † Current address: Department of Plant Molecular Biology and Physiology, Albrecht-von- Haller-Institute for Plant Sciences, University of Göttingen, Julia-Lermontowa-Weg 3, D- 37077 Göttingen, Germany. ¥ Current address: Institut de Pharmacologie et de Biologie Structurale, IPBS, Université de Toulouse, CNRS, UPS, Toulouse, France * To whom correspondence should be addressed: e-mail: nuria.sanchez-coll@cragenomica.es Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB Campus UAB Bellaterrra, 08193 Spain e-mail: david.rengel@ipbs.fr Institut de Pharmacologie et de Biologie Structurale BP 64182 205 route de Narbonne 31077 Toulouse Cedex 04 France Conflict of interest statement: The authors declare no conflict of interest. Abstract Recognition of a pathogen by the plant immune system often triggers a form of regulated cell death traditionally known as the hypersensitive response (HR). This type of cell death occurs precisely at the site of pathogen recognition, and it is restricted to a few cells. Extensive research has shed light into how plant immune receptors are mechanistically activated. However, a central key question remains largely unresolved: how does cell death zonation take place and what are the mechanisms that underpin this phenomenon? Consequently, bona fide transcriptional indicators of HR are lacking, which prevents gaining a deeper insight of its mechanisms before cell death becomes macroscopic and precludes any early or live observation. We addressed this question using the paradigmatic Arabidopsis thaliana– Pseudomonas syringae pathosystem, by performing a spatio-temporally resolved gene expression analysis that compared infected cells that will undergo HR upon pathogen recognition vs by-stander cells that will stay alive and activate immunity. Our data revealed unique and time-dependent differences in the repertoire of differentially expressed genes, expression profiles and biological processes derived from tissue undergoing HR and that of its surroundings. Further, we generated a pipeline based on concatenated pairwise comparisons between time, zone and treatment that enabled us to define 13 robust transcriptional HR markers. Among these genes, the promoter of an uncharacterized AAA-ATPase has been used to obtain a fluorescent reporter transgenic line, which displays a strong spatio-temporally resolved signal specifically in cells that will later undergo pathogen-triggered cell death. In sum, this valuable set of genes can be used to define those cells that are destined to die upon infection with HR-triggering bacteria, opening new avenues for specific and/or high- throughput techniques to study HR processes at a single-cell level. Keywords: Arabidopsis thaliana, Cell Death Indicator, Effector-Triggered Immunity, Hypersensitive Response, Pattern-Triggered Immunity, Plant Immunity, Pseudomonas syringae. 1 INTRODUCTION 2 3 Plants are rich sources of nutrients for pathogens with contrasting lifestyles (Dangl et al., 2013). 4 As opposed to animals, plants do not possess a circulatory system with mobile cells specialized 5 in pathogen defense (Jones and Dangl, 2006). Since their cells are fixed by their cell walls, 6 plants rely on each cell’s autonomous immunity and on systemic signals emanating from 7 infection sites to distal cells to prime the plant for future pathogen encounters (Ausubel, 2005). 8 Moreover, instead of a somatic adaptive immune system that produces antigen receptors on 9 demand, plant cells are equipped with extracellular pattern-recognition receptors (PRRs) and 10 intracellular nucleotide-binding leucine-rich repeat immune receptors (NLRs) that recognize 11 microbe-associated microbial patterns (MAMPs) and pathogen effectors required for 12 virulence, respectively (Couto and Zipfel, 2016). PRR activation brings about a broad defense 13 response named pattern-triggered immunity (PTI), while NLR activation triggers a potentiated 14 and prolonged immune response named effector-triggered immunity (ETI) that reinforce 15 defense outputs observed during PTI (Yuan et al., 2021a; Ngou et al., 2021b). ETI often 16 culminates in a macroscopic localized cell death at the attempted pathogen ingress site known 17 as hypersensitive response (HR)-cell death or immune-related cell death (Olvera-Carrillo et al., 18 2015; Balint-Kurti, 2019; Salguero-Linares and Coll, 2019). 19 20 Regulated cell death has a crucial role in both animals and plant immune responses. Extensive 21 research in the animal field supports the notion that the immune system is highly dependent on 22 cell death for a robust and tightly controlled immune response to occur (Lu et al., 2014; Nagata 23 and Tanaka, 2017) . In plants, our knowledge about the biochemical and genetic pathways 24 regulating cell death, particularly in the context of immunity, is still very limited. As an attempt 25 to shed light into how HR is orchestrated in plants, most efforts have been directed towards 26 understanding how NLRs are mechanistically activated, as well as identifying molecular 27 components upstream or downstream of NLRs that are required for HR to occur (Dangl and 28 Jones, 2019; Wang et al., 2019a; Wang et al., 2019b; Ma et al., 2020; Ngou et al., 2021a) 29 30 Plant NLRs can be broadly classified into TNLs and CNLs based on their domain composition: 31 TNLs contain a Toll/Interleukin-1 Receptor (TIR), whereas CNLs harbor a coiled-coiled 32 domain at their N-terminal end (Jones et al., 2016). Groundbreaking research has shown that 33 in plants, pathogen perception leads to NLR oligomerization, which ultimately will result in 34 cell death and immunity (Wang et al., 2019a; Wang et al., 2019b; Ma et al., 2020; Förderer et 35 al., 2022). Oligomerized forms of CNLs can form pores at the plasma membrane that act as 36 Ca2+-permeable channels (Wang et al., 2019a; Wang et al., 2019b; Jacob et al., 2021). Some 37 TNLs, in turn, can oligomerize upon activation to reconstitute a holoenzyme that triggers cell 38 death by a mechanism that is not fully elucidated but that may involve their NAD+ hydrolase, 39 as well as their 2′,3′-cAMP/cGMP synthetase activities (Ma et al., 2020; Martin et al., 2020; 40 Yu et al., 2021). How oligomerization translates to immune signaling and HR remains to be 41 defined. 42 43 In the context of signaling downstream NLR activation or ETI, large-scale transcriptional 44 studies have highlighted the importance of phytohormone networks for high-amplitude 45 transcriptional reprogramming to mount a fast and efficient response (Mine et al., 2018). 46 Comparisons between host transcriptional responses elicited by PTI and ETI suggest minor 47 qualitative differences in the repertoire of genes differentially expressed (Navarro et al., 2004; 48 Mine et al., 2018). These studies also support the recently evidenced assumption that ETI and 49 PTI share immune signaling components (Pruitt et al., 2021; Yuan et al., 2021a; Ngou et al., 50 2021b). However, a central key question remains unexplored: which early transcriptional 51 signatures differentiate cells that recognize the pathogen and will undergo HR from by-stander 52 cells that will remain alive and will activate defenses to fight the pathogen? In recent literature, 53 a few studies underscore the importance of zonation during HR (Betsuyaku et al., 2018; Giolai 54 et al., 2019; Lukan et al., 2020). At the hormonal level, it has been shown that salicylic acid 55 (SA) plays a major role at pathogen-inoculated spots that will later undergo HR, while the 56 jasmonic acid (JA) signaling pathway is activated in the cells surrounding the central SA-active 57 cells (Dorey et al., 1997; Betsuyaku et al., 2018). Furthermore, precision transcriptomics during 58 the immune response elicited by the potato Ny-1 gene against potato virus Y (PVY) revealed 59 the importance of SA accumulation and genes involved in the generation of reactive oxygen 60 species (ROS) for efficient confinement of macroscopic cell death lesions caused by PVY 61 (Lukan et al., 2020). The cell wall polymer lignin has also been shown to participate in HR 62 zonation, by forming a physical barrier around the infection site upon pathogen recognition 63 that presumably contributes to confine the invading agents and restricts colonization (Lee et 64 al., 2019). A transcriptional meta-analysis of developmental vs HR-cell death in plants could 65 only reveal robust indicators for developmental cell death but not for HR-cell death (Olvera- 66 Carrillo et al., 2015) . We realized that the limitation of previous large-scale transcriptomic 67 analysis lacked the spatial dimension of HR (Lewis et al., 2015; Mine et al., 2018), as dying 68 cells were not compared to by-stander cells, and the focus was not placed on identifying 69 specific cell death markers, but rather bulk-analyzing the ETI response at the inoculated area. 70 71 A systematic gene expression analysis of the zonation of HR overtime would help 72 understanding the process of HR at the molecular level and importantly, would allow defining 73 bona fide transcriptional markers of the process. With this purpose, we generated RNA- 74 sequencing (RNA-seq) data to systematically analyze and compare the transcriptional 75 programs taking place at the zone of inoculation/pathogen recognition that will undergo HR vs 76 the surrounding area that will stay alive and activate immunity. We show unique and time- 77 dependent differences in the repertoire of differentially expressed genes (DEGs) and 78 expression profiles derived from tissue undergoing HR and that of its surrounding tissues. 79 Furthermore, we generated a pipeline based on pairwise comparisons between time, zone and 80 treatment that enabled us to define 13 robust transcriptional HR markers and a fluorescent 81 transgenic reporter line. These valuable set of genes can be used to define those cells that are 82 destined to die upon pathogen recognition before the onset of cell death becomes 83 macroscopically visible, opening new horizons to study the processes therein by live, cell- 84 specific and/or high-throughput techniques. 85 86 RESULTS 87 88 Zonally dissected Arabidopsis transcriptomes upon Pto AvrRpm1 infection reveal unique 89 spatio-temporal gene expression. 90 In our experiments we used the paradigmatic interaction between Arabidopsis thaliana Col-0 91 (hereafter Arabidopsis) and the bacterial pathogen Pseudomonas syringae pathovar tomato 92 (Pto) carrying the effector AvrRpm1 (hereafter Pto AvrRpm1), which triggers restricted HR at 93 the site of inoculation upon recognition by the CNL RPM1 (RESISTANCE TO 94 PSEUDOMONAS SYRINGAE PV MACULICOLA 1) (Mackey et al., 2002). In order to 95 zonally dissect HR and its surrounding, we syringe-infiltrated a limited area (roughly 3-4 mm) 96 at the side edge of Arabidopsis leaves with either a mock solution or Pto AvrRpm1. Collected 97 tissue from this area was designated as the “IN” zone. To ensure proper separation between IN 98 and OUT zones, a buffer zone expanding 1 mm next to the IN area was discarded, and a parallel 99 region expanding 1 to 2 mm towards the vein was designated as “OUT” (Figure 1A). We 100 collected tissue at 0, 1-, 2-, 4- and 6-hours post-inoculation (hpi), extracted RNA and assessed 101 transcript abundance by RNA-seq. Under these conditions, macroscopic cell death started 102 appearing at 4 hpi in the Pto AvrRpm1-inoculated samples, as visualized by trypan blue staining 103 (Figure 1B). As expected, this cell death is concomitant with a dramatic drop in photosynthetic 104 efficiency of photosystem II (Fv/Fm ratio) and electron transport rate (ETR) at the IN area 105 (Figure 1C) (Berger et al., 2007). 106 107 To determine whether the obtained RNA-seq data complied with our working hypothesis of 108 spatio-temporal gene expression regulation we performed a Principal Component Analysis 109 (PCA) (Figure S1A-B). We observed that at the IN area, Pto AvrRpm1-treated samples 110 separated from their mock controls from 2 hpi onwards. At the OUT area, however, only Pto 111 AvrRpm1-treated samples at 4 and 6 hpi separated from mock controls. Overall, the PCA 112 confirms that the biggest changes in gene expression are produced at IN, particularly at 4 and 113 6 hpi, whereas at OUT there is a subtler modulation that is most pronounced at 4 hpi. 114 115 Next, we identified differentially expressed genes (DEGs) between bacteria and mock- 116 inoculated samples (DEGs; false discovery rate (FDR) < 0.05 and |log2FC| > 2), thereby 117 characterizing the transcriptional changes occurring at each tissue area at every time point. We 118 found a total of 5,495 DEGs at the IN zone and 1,785 at the OUT zone (Figure 2A, Table S1). 119 Enrichment of Gene Ontology (GO) terms was examined in every group of DEGs at each 120 specific time point (Figure S2, Table S2). Upregulated genes at the IN area were enriched in 121 immunity- and phytohormone-associated processes (Figure S2A). Immunity-related GO terms 122 associated with PTI and ETI such as “plant-type hypersensitive response” and “pattern 123 recognition receptor signalling pathway” appeared at initial stages of infection (1 and 2 hpi), 124 while at later stages (from 2 hpi onwards) there is enrichment of GO terms associated with 125 more general defense and abiotic stress processes such as “defense response to bacteria” and 126 “response to wounding”, respectively (Figure S2A). Regarding phytohormone-related 127 processes, we observed an enrichment in SA-related GO terms from 1 hpi onwards, confirming 128 the importance of SA at the HR/IN area (Dorey et al., 1997; Zheng et al., 2015). In contrast, 129 GO terms associated with JA were particularly overrepresented at later time points (4 and 6 130 hpi), in accordance with previous findings demonstrating that SA can activate JA signaling 131 through a non-canonical pathway promoting ETI (Liu et al., 2016). GO terms related to other 132 defense/stress-related phytohormones such as ethylene (ET) and abscisic acid (ABA), were 133 also enriched at 4 and 6 hpi (Figure S2A). 134 135 Among downregulated genes at the IN zone, an enrichment in GO terms related to 136 photosynthesis and chloroplast biology occurred at late time points (4 and 6 hpi) (Figure S2B). 137 This correlates with the drop in photosynthetic efficiency shown in Figure 1C, which is part 138 of the defense/yield trade-off to derive resources to immune responses and shut down 139 production of sugars and nutrients, as they might serve as a source for pathogen survival and 140 multiplication (Lu and Yao, 2018). 141 142 Strikingly, at the OUT area we only observed differential expression at late time points (4 and 143 6 hpi), with an overall reduction in the number of DEGs compared to the IN area (Figure 2A). 144 Upregulated genes were enriched in GO terms associated with hormonal regulation, 145 particularly to the JA signaling pathway (Figure S2C). Downregulated genes at the OUT area 146 did not show any enriched GO term, possibly due to the low number of genes. 147 148 To identify genes exclusively upregulated (FDR < 0.05 and |log2FC| > 2) at either the IN or 149 OUT areas we first generated Venn diagrams representing the number of genes modulated at 150 each time point upon infection (Figure S3). This analysis confirmed that upregulation at both 151 IN and OUT mainly occurs at 4 or 6 hpi (Figure S3) and therefore, we selected these two time 152 points to further identify genes that are exclusively upregulated at each tissue area (Figure 2B). 153 Specifically, we found a total of 1,840 genes being upregulated exclusively at IN, 1,117 genes 154 upregulated at both IN and OUT and 221 genes being exclusively upregulated at OUT (Figure 155 2B, Table S3). Among the overrepresented GO terms found in genes exclusive for the IN area 156 were “defense response to bacterium”, “response to molecule of bacterial origin” or “response 157 to salicylic acid”. We also found various GO terms associated with responses to several other 158 stresses such as salt, oxygen-containing compounds, sulfur compounds, heat and hydrogen 159 peroxide (Figure 2C, Table S4), which is not surprising, considering that the tissue is 160 undergoing cell death. In contrast, overrepresented GO terms in genes exclusively upregulated 161 at the OUT area included “regulation of defense response” and, interestingly, “response to 162 wounding” and “response to jasmonic acid” (Figure 2C, Table S4). These JA-related genes 163 follow a very distinct expression pattern, with an early peak at 1 hpi both at the IN and OUT 164 areas, and a second peak at 4 hpi of higher intensity in the OUT zone (Figure S5, Table S4). 165 Although further experimental validation would be required, these data reveal expression 166 patterns of a set of genes that could potentially be used as OUT markers along with previously 167 reported markers such as VSP1 (Chung et al., 2008; Betsuyaku et al., 2018). To better visualize 168 the behavior of the remaining OUT specific genes throughout the course of the infection, we 169 generated heatmaps representing their differential expression at IN and OUT areas (Figure 170 S4). 171 172 Clustering of gene expression profiles reveals distinct expression patterns at the IN and 173 OUT areas over time 174 Next, we set out to determine whether genes at the IN and OUT areas followed specific 175 expression patterns and if particular biological processes were associated to those patterns. For 176 this, we first analyzed gene expression profiles using Fuzzy c-means, a soft partitioning 177 algorithm which offers robust clustering with regards to noise by variation of a fuzzification 178 parameter that limits the contribution of ill-behaved profiles to the clustering process (Olsen et 179 al., 2006; Kumar and Futschik, 2007). Based on this, we could define three and five distinct 180 and non-overlapping clusters for Pto AvrRpm1-treated samples in the IN and OUT areas, 181 respectively (Figure 3, Figure S8 and Table S6-S7). Genes within each cluster were 182 subsequently re-clustered in mock-treated samples, producing two distinct sub-clusters 183 (Figures 3, Figure S8 and Table S6-S7). This procedure provided a more detailed overview 184 of the differences and similarities of trajectories between treatments over time and reflected 185 the well-documented wound response that takes place in mock-treated tissue (Mine et al., 2018; 186 Giolai et al., 2019; Vega-Munoz et al., 2020). 187 188 At the IN area of infection, cluster I exhibited a pattern of upregulation from 0 to 2 hpi and 189 mild downregulation from 2 to 6 hpi (Figure 3A). Genes near to its centroid (MSV > 0.7; see 190 Materials and Methods) are mainly associated with immune-related GO terms (Figure S6A, 191 Table S8). Genes in this cluster followed two distinct trajectories in the mock-treated samples: 192 while mock sub-cluster 1.1 showed a steady increase throughout the experiment, mock sub- 193 cluster 1.2 exhibited a typical wounding immune-related response common with infected 194 samples, peaking at 1 h and rapidly returning to steady-state levels (Savatin et al., 2014). 195 196 Cluster II-IN includes genes with a sharp increase of expression at 4 hpi (Figure 3A). Many 197 of the genes following that trajectory are involved in protein degradation processes (autophagy, 198 protein targeting to the vacuole, proteasome mediated degradation) taking place in response to 199 infection (Figure S6A and Table S8). Sub-clusters from mock-treated samples predominantly 200 followed a similar steady trajectory throughout the experiment, which points to an infection- 201 specific effect of upregulation on protein turnover due to infection at the IN area (Figure 3A, 202 Figure S7A and Table S10). 203 204 Cluster III-IN exhibits an expression pattern of steady downregulation from 0 to 4 hpi, followed 205 by a slight recovery of expression from 4 to 6 hpi (Figure 3A). This cluster includes mostly 206 genes belonging to GO terms related to photosynthesis (Figure S6A and Table S8). In this 207 case, mock-treated samples sub-cluster into two distinct patterns of expression: sub-cluster 3.1 208 follows a similar pattern as infected samples, while sub-cluster 3.2 shows a transient decrease 209 of expression at 1 h followed by a recovery phase from 2 to 6 hpi (Figure 3A). Our data show 210 that only certain components of the photosynthetic machinery are specifically affected by the 211 pathogen treatment (Figure S6A, Figure S7A and Table S10). 212 213 At the OUT area of infection, cluster I includes genes that display a sharp peak of expression 214 at 4 hpi (Figure 3B). From this cluster, genes near the centroid belong to GO terms associated 215 with metabolism, hormonal regulation, and wounding response, among others (Figure S6B 216 and Table S9). Interestingly, JA- and SA-responsive genes, which are known to act 217 antagonistically and cooperatively during ETI (Liu et al., 2016; Betsuyaku et al., 2018), seem 218 to be highly enriched in the OUT area. Genes comprising the mock-derived sub-clusters follow 219 a similar trend of steady expression throughout the time course of the experiment, suggesting 220 that the peak of high expression is a specific response to the bacterial infection in the 221 surrounding area (Figure 3B, Figure S7B and Table S11). 222 223 Cluster II-OUT in Pto AvrRpm1-treated samples follows an expression pattern with two sharp 224 upregulation peaks at 1 and 4 hpi (Figure 3B). These trajectories are followed by genes 225 associated with JA-related processes and wounding, and is a very specific pattern exclusively 226 found at the OUT zone (Figure 3, Figure S6B and Table S9). The early peak at 1 hpi shared 227 between mock and infected samples could account for a wounding response elicited early at 228 the area surrounding the syringe-infiltrated area, whilethe peak at 4 hpi appears as a late 229 response that occurs specifically at the tissues surrounding the pathogen inoculation area 230 (Figure 3, Figure S6 and Table S11). 231 232 In cluster III-OUT, the trajectory of genes from Pto AvrRpm1-treated samples does not 233 remarkably differ from mock treatment (Figure 3B). Genes that comprise this cluster mainly 234 fall into GO terms associated with the photosynthetic machinery (Figure S6B and Table S9). 235 These data indicate that photosynthesis at the OUT area of infection does not seem to be altered 236 by pathogen infection as opposed to the IN area (Figure 3B and Figure S6-S7) correlating 237 with zonal photosynthesis efficiency values shown in Figure 1C and as previously reported 238 (Berger et al., 2007). 239 240 Novel zonal HR transcriptional indicators can be elucidated from pairwise comparisons 241 between time, treatment and area 242 In order to identify robust HR markers that are exclusively upregulated at the site of cell death 243 (IN area) we conducted a pipeline of differential expression analysis that consisted of 244 concatenated pairwise comparisons considering the three variables in our experimental design: 245 time, treatment and area (Figure 4A). Since the highest degree of differential expression 246 between treatments took place at 4 and 6 hpi (Figure 2A), we carried out the comparisons at 247 these two time points independently. Firstly, we focused on the time variable and selected 248 genes that were confidently upregulated at the IN area of Pto AvrRpm1-infected plants, either 249 at 4 and/or 6 hpi, compared to 0 hpi (1st filter: FDR < 0.05 and log2FC>2). Secondly, we 250 removed genes also upregulated at 4 and/or 6 hpi at the IN area in mock controls (2nd filter: 251 FDR < 0.05 and log2FC>2). Since we aimed to find genes only upregulated at the IN/cell death 252 area, next, we removed the genes that were upregulated by bacterial inoculation at the OUT 253 area at least to half of the levels than in the IN zone (3rd filter: FDR <0.05 and log2FC<1). 254 Finally, from the genes that met those three criteria, we kept those that were differentially 255 upregulated at IN compared to the OUT area in Pto AvrRpm1-infected plants (4th filter: FDR 256 < 0.05 and log2FC > 2) (Figure 4A). 257 258 A total of 32 genes passed all 4 filters, constituting a set of potential HR indicators (Figure 259 S9). From these, 24 were extracted from the 4 hpi dataset and 11 from the 6 hpi dataset and 3 260 from both time points (Figure S9). Of note, due to the stringency of the filters none of these 261 genes passed all the filters at 1 or 2 hpi, although 7 of them were upregulated after infection at 262 the IN zone at these early time points (M5, M6, M7, M8, M9, M11 and M12). The expression 263 profiles of this putative HR indicators can be visualized as DESeq2 pseudo-counts as a function 264 of time at both areas of infection in Figure S10. The expression patterns of these 32 genes at 265 0 and 4/6 hpi were validated by real time quantitative PCR (RT-qPCR) using newly obtained 266 biological samples (Figure S11). To ensure that the potential markers were exclusively 267 upregulated as part of the HR response triggered by effector-mediated bacterial recognition 268 and not as part of the defense responses triggered by disease-causing bacteria, we also included 269 samples inoculated with Pto DC3000 EV (Pto EV), a strain that causes disease but does not 270 trigger HR in Arabidopsis Col-0. Among the 32 genes tested, a total of 13 (10 of them at 4 hpi, 271 4 at 6 hpi, with one at both time points), behaved as bona fide HR indicators (Figure 4B-C), 272 showing a distinctive upregulation specifically triggered at the IN area by an HR-causing 273 bacterium. 274 275 The AAA-ATPase At5g17760 promoter specifically drives expression of GFP to the IN 276 area of infection, constituting a robust transcriptional live marker of HR 277 In order to generate much needed tools to extend our understanding of how HR unfolds at the 278 infection site and its surrounding tissue, we generated stable transgenic Arabidopsis plants 279 expressing green fluorescent protein (3xGFP) under the control of the promoters of each of the 280 13 identified putative HR marker genes. A nuclear localization signal (NLS) was fused to GFP 281 to concentrate the signal in the nucleus and facilitate detection, which enabled us to distinguish 282 promoter-driven fluorescence from the auto-fluorescence derived from HR (Bennett et al., 283 1996). 284 285 We focused our analysis on plants expressing pAT5G17760:NLS-3xGFP (corresponding to 286 M13), as they showed high, cell-specific, robust and clear GFP signal in the nuclei of the leaf 287 regions infected with Pto AvrRpm1 (Figure 5B, Figure S12). In several independent transgenic 288 lines, activation of pAT5G17760 was limited to the syringe-infiltrated area and could not be 289 detected in the surrounding tissues (Figure S13). In all pAT5G17760:NLS-3xGFP marker lines 290 the GFP signal appeared concomitantly with cell death, as shown by trypan blue staining 291 (Figure 5B, Figure S12 and S13). A clear GFP signal was not detected in all other marker 292 lines tested. 293 294 In addition to Pto AvrRpm1, we also analyzed the response of pAT5G17760:NLS-3xGFP plants 295 to Pto expressing AvrRpt2 (Pto AvrRpt2), which induces HR in Col-0 plants via the CNL 296 RESISTANT TO P. SYRINGAE 2 (RPS2) (Mackey et al., 2003) and to Pto expressing 297 AvrRps4 (Pto AvrRps4), where HR is mediated by the TNL-pair RPS4/RRS1 and requires 298 helper NLRs (Gassmann et al., 1999; Narusaka et al., 2009). The same pattern was observed 299 after infiltration with Pto AvrRpt2 or Pto AvrRps4 (Figure 5B), which indicates that 300 pAT5G17760 robustly responds to pathogen-mediated activation of different classes of NLR 301 receptors. As controls, we included mock, Pto EV and a non-pathogenic mutant strain secreting 302 no effectors (Pto hrcC-) (Alfano et al., 2000). Importantly, infiltration with the mock solution 303 or with non-HR causing bacterial strains did not activate pAT5G17760. It is worth noting that 304 for microscopy imaging experiments we used a lower bacterial inoculum (O.D600 0.01) to 305 mimic more natural infection conditions and to delay the onset of HR and tissue collapse 306 (Figure 5A), which was necessary for microscopic detection of GFP. At higher inoculum, rapid 307 accumulation of phenolic compounds at the site of infection results in extremely high auto- 308 fluorescence levels that hamper imaging. 309 310 Since pathogens with contrasting lifestyles can trigger HR or HR-like cell death in plants, we 311 tested whether this reporter line can be employed in a broader sense. For this, we infected adult 312 Arabidopsis leaves by drop inoculation with Botrytis cinerea, a necrotrophic pathogen that kills 313 plant tissue prior to feeding, using a range of toxic molecules (Muckenschnabel et al., 2002). 314 At 3 dpi, we observed GFP expression in the nuclei of cells at the region inoculated with the 315 pathogen as opposed to mock-inoculated plants (Figure S14). Together, our observations 316 indicate that pAT5G17760 activity is spatially regulated and confined to the area undergoing 317 HR elicited by hemibiotrophic (P. syringae) and necrotrophic (B. cinerea) pathogens. Thus, 318 the transgenic reporter line pAT5G17760:NLS-3xGFP constitutes a very useful tool to monitor 319 this process in planta. 320 321 The AT5G17760 gene encodes a putative AAA ATPase of unknown function. A knock-out 322 mutant of this gene did not show any alteration in HR or pathogen growth restriction compared 323 to wild-type plants (Figure S14). The lack of phenotype could be due to functional 324 redundancy/compensation, a very common masking phenomenon in plants. 325 326 HR markers and particularly At5g1776 are highly upregulated in other RNA-seq data 327 sets from plants undergoing ETI and autoimmunity 328 We looked at the behavior of At5g17760 and the rest of marker genes in already published 329 RNA-seq data sets from either plants undergoing ETI or autoimmune mutant plants displaying 330 constitutive defense responses and runaway cell death (Figure S16) (Mine et al., 2018; Yang 331 et al., 2020; Barragan et al., 2020; Chantarachot et al., 2020). Fold changes from marker genes 332 with significant p values (FDR < 0.05) in these data sets were plotted as heatmaps to reveal 333 their level of upregulation (Figure S16). As expected, most gene markers are significantly 334 (FDR < 0.05) upregulated during ETI triggered by Pto AvrRpm1 and Pto AvrRpt2 at 4, 6 and 335 9 hpi in Mine et al., (2018) (Figure S16A). Interestingly, At5g17760 is the highest upregulated 336 gene in hos15-4 and rh6812 mutant plants undergoing autoimmunity (Yang et al., 2020; 337 Chantarachot et al., 2020). Likewise, upregulated genes from data sets of incompatible 338 Arabidopsis F1 hybrids (Cdm-0 x TueScha-9) exhibiting autoimmunity comprised most HR 339 markers found in this study with At5g17760 being the highest upregulated gene (Barragan et 340 al., 2020)(Figure S16B). 341 342 DISCUSSION 343 Zonation of HR in plants is underscored by distinct gene expression patterns and 344 processes in dying vs by-stander cells 345 In plants, pathogen recognition via intracellular NLR receptors often results in an HR reaction 346 that helps preventing pathogen proliferation (Pitsili et al., 2020). This is a highly zonal response 347 that takes place at the site of infection, whereby dying cells send signals to the surrounding 348 tissues to activate defenses and block pathogen invasion. Traditionally, the plant immune 349 system was considered strictly two-branched, with PTI being elicited by recognition of 350 conserved pathogen patterns via cell surface receptors, and ETI recognizing pathogen effector 351 proteins secreted into the plant cell via intracellular NLR receptors (Jones and Dangl, 2006). 352 Over the last decades, many efforts have been directed towards understanding the 353 transcriptional reprogramming elicited during PTI and ETI (Tao et al., 2003; Lewis et al., 2015; 354 Bozso et al., 2016; Mine et al., 2018; Duan et al., 2020). One of the major conclusions drawn 355 from these studies is that whilst the repertoire of differentially expressed genes in the host is 356 largely similar, ETI leads to a faster and more robust transcriptional response than PTI (Tao et 357 al., 2003; Mine et al., 2018; Yuan et al., 2021a; Yuan et al., 2021b; Ngou et al., 2021b). These 358 findings, together with emerging evidence showing additional levels of synergy and crosstalk 359 between PTI and ETI has somewhat blurred the traditional PTI-ETI dichotomy (Dongus and 360 Parker, 2021; Ngou et al., 2021a; Pruitt et al., 2021). However, despite the large amount of 361 time-resolved transcriptomic data produced (Tao et al., 2003; Lewis et al., 2015; Hillmer et al., 362 2017; Mine et al., 2018), the spatial consideration of HR upon ETI activation has been partly 363 overlooked, with only few studies pointing to its importance in regulating the process (Dorey 364 et al., 1997; Betsuyaku et al., 2018; Giolai et al., 2019; Lukan et al., 2020). It remains unclear 365 whether and to what extent transcriptional reprogramming takes place at the vicinity of cell 366 death compared to that occurring at the infected area upon bacterial infection. 367 368 Our experimental design (Figure 1A) considered the spatio-temporal angle of plant HR to gain 369 a better understanding of how this process is restricted to a few cells upon pathogen recognition 370 and to define robust markers of the dying area over time. This is particularly important since 371 in plants, cell death characterization has largely relied on biochemical and morphological 372 hallmarks most of which are post-mortem and which in most cases do not provide unequivocal 373 criteria (van Doorn, 2011; van Doorn et al., 2011). We currently lack a set of genes that can be 374 employed as gene indicators of cell death triggered by pathogens. In silico comparisons of 375 transcriptome profiles at different developmental stages and upon environmental stresses 376 leading to cell death, enabled identification of cell death indicators of developmentally 377 regulated programmed cell death that can be used to detect or even isolate cells that are ready 378 to die (Olvera-Carrillo et al., 2015). The same approach did not lead to identification of reliable 379 HR markers, partly because the available datasets were not obtained on zonally resolved 380 samples (Olvera-Carrillo et al., 2015). 381 382 Here, differential expression analysis and clustering of genes based on expression profiles over 383 time enabled us to infer biological processes taking place at each tissue area (IN/OUT) upon 384 bacterial infection, giving us hints on how HR can be spatially restricted. At the IN area, genes 385 involved in a local immune response to ETI-triggering bacteria are greatly induced from 1 hpi 386 onwards (cluster I) (Figure 2A and Figure S2A, Figure 3A). Tissue from the IN area, also 387 contains a set of genes that show a peak of upregulation from 2 to 4 hpi (cluster II), involved 388 in diverse biological processes ranging from regulation of immunity, responses to JA and SA 389 or protein turnover (Figure S6A). It is now well established that proteasome activity is strongly 390 induced during bacterial infection and that certain subunits of the proteasome are required for 391 efficient fine-tuning of immune responses in plants (Misas-Villamil et al., 2013; Ustun et al., 392 2016; Ustun et al., 2018). Finally, we identified a strong transcriptional repression of 393 photosynthetic genes at 4 hpi at the IN area (cluster III) (Figure 2B, Figure S2B, Figure 3A 394 and Figure S6A), in accordance with the previously established notion that infection results 395 in a global downregulation of genes associated with the photosynthetic machinery (Bilgin et 396 al., 2010). This specific decrease in photosynthesis is particularly interesting in light of recent 397 reports of the interplay between bacterial effectors and the chloroplast, whereby certain 398 effectors can suppress chloroplast functions and in turn, chloroplasts can adopt immune 399 functions to fight off pathogens (Kachroo et al., 2021; Littlejohn et al., 2021; Savage et al., 400 2021). 401 402 Our results also show that transcriptional reprogramming in host cells surrounding the infection 403 area (OUT area) is less extensive with a lower number of differentially expressed genes than 404 at the IN area, and starts later mostly from 4 hpi onwards (Figure 2A). Remarkably, 405 photosynthesis is not significantly affected at the OUT area, corroborating our in vivo 406 measurements (Figure 1C) and previous findings (Bilgin et al., 2010). A relatively functional 407 photosynthetic machinery may be key to maintain effective defense mechanisms and prevent 408 these cells from dying as their neighbors. This finding might have been masked in previous 409 transcriptional studies that have not taken into account the zonal nature of HR, and reveals that 410 the defense-growth trade-off may also have a marked spatial component that needs to be taken 411 into account in future research. Besides photosynthesis, the OUT zone was characterized by a 412 marked upregulation of wound/JA-related genes at 4 hpi (Figure 2C, Figure 3B and Figure 413 S2C). This response can also be observed at the IN zone but the level of upregulation at the 414 OUT zone is remarkably higher (Figure S4), indicating an amplification in JA signaling at the 415 cells surrounding the death zone. In addition, some of the JA-related genes are among those 416 genes exclusively upregulated at OUT at 4/6 hpi, which indicates that they could potentially 417 be used as zonal markers of the surrounding area (Figure 2B-C and Figure S5). In vivo 418 imaging of marker gene promoter activities of SA and JA signaling during ETI discerned two 419 spatially distinct domains around the infection site, where JA signaling is thought to be 420 important for regulating over-activation of SA signaling (Betsuyaku et al., 2018). Future 421 studies that include mutants deficient in JA could provide mechanistic insights into how JA 422 signaling contributes to the confinement of plant HR. Our analysis also shows that some SA- 423 signaling genes are among the upregulated IN-specific genes at late time points (Figure 2B-C 424 and Table S4). Although originally considered antagonistic hormones required for immunity 425 against pathogens with contrasting lifestyles (Spoel et al., 2007), the interplay and synergism 426 of these two phytohormones is now well established during ETI (Liu et al., 2016). 427 428 Zonally resolved transcriptomic analysis of HR in plants allows for the identification of 429 robust biomarkers of the process 430 431 Robust biomarkers are essential to gain mechanistic knowledge of cell- or tissue-specific 432 processes. In mammals, the extensive mechanistic knowledge of molecular constituents 433 underlying regulated cell death has enabled the use of biomarkers for detection of tumor cells 434 or aberrant cell death processes in cancer patients (Abu-Qare and Abou-Donia, 2001; Ward et 435 al., 2008). The field of HR in plants is gaining momentum thanks to recent major discoveries 436 that in one hand are leading to a redefinition of the PTI-ETI relationship and on the other, have 437 provided mechanistic insight into how NLRs become activated and form supramolecular 438 complexes that mediate cell death (Wang et al., 2019a; Wang et al., 2019b; Martin et al., 2020; 439 Bi et al., 2021; Jacob et al., 2021; Pruitt et al., 2021; Tian et al., 2021; Yuan et al., 2021a; Ngou 440 et al., 2021b; Förderer et al., 2022). However, amidst this exciting scenario, the conceptual 441 framework of HR zonation is scarcely defined and will be key to understand its execution, 442 spatial restriction mechanisms and define bona fide indicators of the process. 443 444 One of the main goals of our analysis was to define new markers of HR. We made use of the 445 RNA-seq data generated from IN and OUT areas to pinpoint gene indicators of HR that can be 446 used either as transcriptional markers or gene promoter markers for in planta detection of cells 447 destined to die using live imaging. Applying stringent filters to our dataset we identified 13 448 genes that can be used as unequivocal transcriptional markers of zonally restricted cells that 449 have activated a death program in response to pathogen perception via NLR activation (Figure 450 4C). 451 452 This marker set includes genes involved or putatively involved in various processes such as 453 ion transport across the plasma membrane (M1), cell detoxification (M2, M3), lipid 454 metabolism (M5, M6), cell wall remodeling (M7, M8, M9), protein degradation (M10), 455 glycolysis (M11, M12), whereas one of these genes remains largely uncharacterized (M13) and 456 encodes an AAA+ ATPase of unknown function. Interestingly, all these predicted functions 457 are consistent with processes expectedly taking place on cells destined to die or that have 458 started dying, although the function of most of these genes remains to be fully determined. This 459 set of genes provide a glimpse of transcriptional regulation of HR at the site of infection, the 460 tip of the iceberg of the multi-level regulation of the process. For example, the fact that several 461 genes are involved in cell wall remodeling highlights the importance of processes taking place 462 in this extracellular compartment. In line with this, an increase of lignification at the edge of 463 cells undergoing HR was shown in the past and provided a clear picture of the zonal nature of 464 this process (Lee et al., 2019). Interestingly, our transcriptome data clearly shows that many 465 lignin biosynthetic genes are strongly and specifically upregulated at the IN zone at certain 466 time points (Figure S15). How this cell wall lignification is regulated upon pathogen 467 perception remains to be clarified and will be an interesting topic of research in the future. 468 469 Our data also reinforces the idea that the proteases involved in degradation of cell components 470 during HR are not particularly regulated at the transcriptional level. We observe specific 471 upregulation of degradative processes at the IN zone such as autophagy, vacuolar degradation, 472 and proteasome-mediated processes and in fact, one of the marker genes is a proteasome 473 subunit (Figure 3 and Figure 4B). However, we did not find any protease specifically 474 upregulated at the IN zone, nor did any of them pass the filters that constitute a marker gene in 475 our study. 476 477 In parallel, the changes observed in marker genes involved either in ion transport across the 478 plasma membrane or cell detoxification may be somewhat related with the predicted formation 479 of a pore at the plasma membrane by pathogen-mediated activation of certain NLRs (Bi et al., 480 2021; Jacob et al., 2021; Förderer et al., 2022). Although crucial pieces of this mechanism have 481 been unveiled, knowledge is still scattered and we lack a more integrated picture that combines 482 NLR activation with downstream processes, including cell death execution. Interestingly, 7 out 483 of the 32 gene markers (M5, M6, M7, M8, M9, M11 and M12) that pass our filters exhibit an 484 early upregulation at 2 hpi compared to mock controls according to the RNA-seq data (Figure 485 S10). Although these genes did not pass the stringent 4-tier filtering applied (Figure 4A) at 2 486 hpi, the expression profiles of these genes could be compatible with their potential use as earlier 487 markers of HR at the IN area. 488 In sum, our data provides a snapshot of how infected cells respond to pathogen recognition at 489 the transcriptional level, compared to their neighbors, that are not directly exposed to the 490 pathogen but respond to it. Importantly, this analysis has revealed a set of genes that are 491 specifically upregulated at the IN zone and constitute robust markers of HR, opening new paths 492 to deepen our knowledge on the process. 493 494 Importantly, we present an Arabidopsis HR reporter line stably expressing GFP under the 495 control of the AAA+ ATPase At5g17760 (M13), which shows extremely clear and strong 496 expression exclusively at the inoculated area where pathogen recognition takes place via ETI, 497 before the onset of cell death becomes apparent (Figure 5A-B and Figure S12). The other 498 genes (M1-M12) constituted very clear qPCR markers but GFP promoter fusions did not result 499 in clear GFP expression. This can be attributed to the limitations from defining an active 500 promoter sequence. 501 502 Interestingly, expression of the marker pAt5g17760:NLS-3xGFP is similarly regulated by 503 different classes of NLRs (CNLs and TNLs) revealing conservation of the process (Figure 504 5A). Moreover, the marker is also induced zonally by necrotrophic pathogens such as B. 505 cinerea that cause an HR-like phenotype (Figure 5C). Thus, this transgenic line constitutes a 506 robust in planta biomarker of HR triggered by activation of different NLRs upon infection with 507 pathogens with contrasting lifestyles. 508 509 Future in-depth analysis of all HR marker genes identified in this work, including 510 combinatorial genetics, will contribute to a better understanding of HR. This set of genes 511 constitutes an invaluable tool to zonally discriminate cells undergoing pathogen-triggered cell 512 death and mechanistically dissect this process. Of particular interest will be to sort GFP- 513 expressing cells of the pAt5g17760:NLS-3xGFP transgenic line upon infection and adapt high- 514 throughput cell death monitoring equipment used so far for animal cell death to describe and 515 quantify the features and regulatory networks that define HR in plants at a single-cell level. 516 517 MATERIALS AND METHODS 518 Plant and bacteria materials and growth 519 The Arabidopsis thaliana accession Col-0 was used for all experiments carried out in this study 520 expect for electrolyte leakage. For electrolyte leakage, Col-0, rpm1-3 (Grant et al., 1995) 521 mutant of the NLR RPM1 and at5g17760 mutant (GABI-KAT line 592F04_1) which carries 522 T-DNA insertion in exon two, were used. Primers used for identifying the T-DNA mutant and 523 for corroboration null expression by RT-qPCR are listed in Table S12. 524 525 Seeds were sown on ½ Murashige and Skoog (MS) media supplemented with 1% sucrose and 526 stratified at 4ºC for two days. Plants were grown in a controlled chamber with a photoperiod 527 of 9 h light and 15 h dark with white fluorescent lamps under 65% relative humidity. Seeds 528 were germinated on plates and grown for 10-7 days, then individually transplanted to Jiffy 529 pellets and grown for 3 additional weeks. 530 531 Pseudomonas syringae pathovar tomato (Pto) strains Pto AvrRpm1, Pto AvrRpt2, Pto AvrRps4, 532 Pto hrpC- and Pto empty vector pVSP61 (EV) were grown on selective King’s B (KB) medium 533 plates for 48 h at 28 °C. Bacteria was then resuspended in 10 mM MgCl2 and the OD600 adjusted 534 to the appropriate inoculum. 535 536 Bacterial inoculation and RNA-seq data collection. 537 Bacteria were resuspended and the concentration was adjusted at 5*107 colony-forming units 538 or to an optical density measured at a wavelength of 600 nm (OD600) of 0.05. Fully expanded 539 7th or 8th rosette leaves were used for infiltration with either a mock solution (10 mM MgCl 2) 540 or Pto AvrRpm1. We syringe-infiltrated an area of roughly 3-4 mm at the side edge of leaves. 541 Upon infiltration, the edge of the infiltrated area was underlined using India ink, and the total 542 area infiltrated designated as “IN”. A 1 mm buffer zone next to the IN area was discarded and 543 used as a reference to properly separate between the IN and the OUT zone, that expanded 1-2 544 mm towards the vein. Leaf tissue was separately collected from the IN and OUT area of 545 infiltration at 5 different time points: 0, 1, 2, 4 and 6 hours by making use of a sterile scalpel. 546 Leaf tissue was stored in 2 mL Eppendorf tubes and snapped-frozen in liquid nitrogen until the 547 time of RNA extraction. Each sample collected consisted of tissue from six leaves derived from 548 three different plants. For generation of three biological replicates from each condition (area, 549 treatment and time), three independent experiments were performed. total sum of 60 samples 550 -2 treatments (mock/infected), 5 time points (0, 1, 2, 4 and 6 hpi), 2 areas (IN/OUT) and 3 551 biological replicates- were used for RNA-sequencing. 552 553 For RNA library preparation, 1 μg of RNA from each sample was isolated using the 554 NucleoSpin® RNA isolation kit (Macherey-Nagel, Hoerdt Cedex, France) following the 555 manufacturer’s instructions. RNAseq was performed at the GeT-PlaGe core facility, INRA 556 Toulouse. RNA-seq libraries have been prepared according to Illumina’s protocols using the 557 Illumina TruSeq Stranded mRNA sample prep kit to analyze mRNA. Briefly, mRNA was 558 selected using poly-T beads. Then, RNA was fragmented to generate double stranded cDNA 559 and adaptors were ligated to be sequenced. 11 cycles of PCR were applied to amplify libraries. 560 Library quality was assessed using a Fragment Analyzer and libraries were quantified by qPCR 561 using the Kapa Library Quantification Kit (Kapa Biosystems, Inc, Wilmington, MA, USA). 562 RNA-seq experiments have been performed on an Illumina HiSeq3000 using a paired-end read 563 length of 2x150 bp with the Illumina HiSeq3000 sequencing kits. 564 565 Read mapping and differential expression analysis 566 “FastQC” and “TrimGalore!” software was used for raw Illumina reads quality control analysis 567 and trimming of reads containing adaptor- or vector-derived sequences, respectively 568 (Babraham Bioinformatics - FastQC A Quality Control tool for High Throughput Sequence 569 Data, 2021). rRNA was detected and removed using “SortMeRNA 2.1b” software (Kopylova 570 et al., 2012). Cleaned reads together with the transcriptome of Arabidopsis thaliana (as of 30 571 August 2018), including ncRNA, were used to quantify gene expression at transcript level 572 using the software “Salmon v0.11.3” (Patro et al., 2017). Raw counts aggregated by gene were 573 obtained using “tximport v1.14.2” and the result was used as input to “DESeq2” v1.26.0 (Love 574 et al., 2014; Soneson et al., 2015) to perform differential expression analysis. Then, genes 575 adding up to less than 10 counts across all 60 samples were removed. The pre-filtered DESeq2 576 object contained 32,865 rows that turned to 23,986 after filtering. Counts normalized for 577 sample size and regularized-logarithm transformed were used to produce PCAs. 578 579 Raw counts together with sample size information were used as input for DESeq2’s differential 580 expression analysis. Simple pairwise comparisons based on a single factor were performed 581 using DESeq2’s “result” function while time course differential expression results were 582 obtained using a likelihood ratio test as previously described (Love et al., 2015). Genes with 583 FDR below 0.05 and |log2FC| higher than 2 were considered as differentially expressed. FDR 584 was calculated according to the Benjamini and Hochberg’s (BH) method (Benjamini and 585 Hochberg, 1995). 586 587 Gene clustering 588 Gene clustering was performed using Mfuzz v2.46.0 package under the R environment (Kumar 589 and Futschik, 2007; RStudio | Open source & professional software for data science teams, 590 2021) which is based on fuzzy c-means clustering algorithms. IN and OUT samples were 591 independently analyzed. After time course differential expression analysis using DESeq2, only 592 genes with an FDR <0.05 in the likelihood ratio test were selected for clustering. 593 594 The optimal number of non-overlapping clusters with a correlation value below 0.85 was 3 and 595 6 for Pto AvrRpm1-treated samples at the IN and OUT areas of infection, respectively. 596 Subsequently, two highly redundant clusters were merged for OUT samples, yielding 5 final 597 clusters. Genes that integrated each cluster derived from Pto AvrRpm1-treated samples were 598 re-clustered for mock-treated samples in order to inspect the differences and similarities of 599 trajectories between treatments over time. Between two and four mock-based sub-clusters were 600 obtained for every infected-cluster. To avoid overlap, we reduced the number of sub-clusters 601 to two in mock-treated samples. Each gene belonging to a cluster returned an associated 602 membership score value (MSV) that ranged from 0 to 1 depending on how well it fitted the 603 expression profile dictated by the overall genes comprising the cluster. Genes that integrate 604 each cluster in Figure 3 area found in Tables S6-S7. 605 606 Enriched Gene Ontology analysis. 607 The set of genes that belonged to expression profile clusters or that exhibited differential 608 expression were input into TAIR for Gene Ontology enrichment analysis for biological 609 processes, which uses the PANTHER Classification system that contains up to date GO 610 annotation data for Arabidopsis (Berardini et al., 2004). The most specific term belonging to a 611 particular family of GO terms was always selected for plotting. Only those GO terms exhibiting 612 an FDR < 0.05 after Bonferroni Correction for multiple testing and a fold enrichment above 2 613 were selected for representation in dot plots. 614 615 Identification of HR indicators 616 For identification of HR indicators, we concatenated four pairwise comparisons using DESeq2, 617 in which we set different thresholds of log2FC, while keeping a stringent cut-off of FDR <0.05 618 throughout all comparisons. Briefly, we firstly selected genes that were upregulated (log2FC 619 > 2) after Pto AvrRpm1 infection at 4 or 6 hpi vs 0 hpi. From the genes that complied with this 620 first filter, we selected those that were specifically upregulated in Pto AvrRpm1-infected vs 621 mock-inoculated samples at 4 or 6 hpi (log2FC >2). From the genes that passed these two filters 622 we kept those with a log2FC <1 at the OUT area in Pto AvrRpm1-infected vs mock-inoculated 623 samples at 4 or 6 hpi. Since genes with log2FC near 0 do not usually have a low FDR, we kept 624 our stringent FDR threshold while setting the log2FC threshold below 1 in order to capture 625 with statistical confidence downregulated and only mildly upregulated genes at this tissue area. 626 Finally, from the genes that met those three criteria, we kept those that were differentially 627 upregulated at the IN area compared to the OUT area in Pto AvrRpm1-infected plants. 628 629 Validation of gene expression by real time quantitative PCR. 630 The same experimental setup used for RNA-seq data generation was followed for experimental 631 validation by RT-qPCR including infections with Pto AvrRpt2, Pto AvrRps4, Pto hrpC- and 632 Pto EV. Briefly, tissue was snap frozen and RNA isolated with the Maxwell® RSC Plant RNA 633 kit (Promega). 1 µg of RNA was reverse transcribed into cDNA with the High-Capacity cDNA 634 Reverse Transcription Kit with RNase inhibitor (Applied Biosystems TM). RT-qPCRs were 635 performed with LightCycler® SYBRgreen I master (Roche) in a LightCycler® 480 System 636 (Roche). Data was analyzed using the ΔΔCT method and represented as fold enrichment of the 637 time point tested (4 or 6 hpi) relative to 0 hpi. Primers for RT-qPCR used in this study are listed 638 in Table S12 along with primer concentrations. RT-qPCR results in numeric format along with 639 Cp values of Targets and Cp values of reference housekeeping gene are listed in Table S13. 640 641 Cell death analysis 642 Trypan blue staining of Arabidopsis leaves was performed by collecting whole leaves in 50 ml 643 tubes (each leaf in a separate tube) at the specified time-points after treatment and covered with 644 a 1:3 dilution of the stain. Tubes were incubated in previously boiled water for 15 min, and 645 then cleared overnight with chloral hydrate on an orbital shaker. After removal of staining 646 solution, leaves were covered in a 50% glycerol solution and photographed using a Leica DM6 647 microscope. 648 649 Electrolyte leakage 650 Whole leaves from four to five-weeks-old Arabidopsis Col-0, rpm1-3 or at5g17760 (GABI- 651 KAT: 592F04) grown in short-day with a photoperiod of 9h light and 15h dark, were infiltrated 652 with Pto AvrRpm1 at a wavelength of 600 nm (OD600) of 0.05 using a 1-ml needleless syringe. 653 Leaf discs were dried and subsequently collected with a 0.8-cm-diameter cork borer from 654 infiltrated leaves. Discs were washed in deionized water for 1 h before being floated on 2 ml 655 deionized water. Electrolyte leakage was measured as water conductivity with a pocket water 656 quality meter (LAQUAtwin-EC-11; Horiba, Kioto, Japan) at the indicated time points. 657 Bacterial growth assay 658 Whole leaves from four to five-weeks-old Arabidopsis Col-0, rpm1-3 or at5g17760 (GABI- 659 KAT: 592F04) grown in short-day conditions (9h light and 15h dark) were infiltrated with Pto 660 AvrRpm1 at a wavelength of 600 nm (OD600) of 0.001 using a 1-ml needleless syringe. Two 661 leaf discs from two different leaves were collected using a 6 mm-diameter cork borer (disc area 662 0.282 cm2). Samples at day 0 and day 3 post-infection were grounded in 10 mM MgCl2 and 663 serially diluted 5, 50, 500, 5,000 and 50,000 times on a 96-well plate. Subsequently, dilutions 664 were spotted (10 μl per spot) on KB medium with antibiotics. The number of colony forming 665 units (CFUs) per drop was calculated and bacterial growth represented as log10 CFU per cm2 666 of tissue. 667 668 Chlorophyll fluorescence imaging 669 An IMAGING-PAM (Pulse-Amplitude-Modulated) M-Series Chlorophyll Fluorometer 670 system (Heinz Walz, Effeltrich, Germany) was used to investigate spatio-temporal changes in 671 photosynthetic parameters at the IN and OUT areas of infection (Schreiber, 2004). Plants were 672 kept in the dark for 30 minutes before measurement. Plants were exposed to 2 Hz frequency 673 measuring light pulses for Fo (minimum fluorescence in the dark-adapted state) determination. 674 Saturating pulses (800 ms) of white light (2400 mmol photons.m-2 s-1) were applied for Fm 675 (maximum fluorescence in the dark-adapted state) determination. The photosynthetic 676 efficiency or maximum quantum yield of PSII photochemistry (Fv/Fm) was determined as 677 (Fm-Fo)/Fm. The relative PSII electron transport rate (ETR) was calculated by performing a 678 kinetic analysis for 10 minutes with 60 second pulses (Schreiber et al., 2012). Areas of interest 679 (AOI) included IN and OUT to evaluate spatial heterogeneity. The measurements were taken 680 after 0, 1, 2, 4 and 6 hpi. Results are shown from 6 different AOI. 681 682 Generation of transgenic promoter reporter lines 683 Regions of approximately -2.5 kb upstream of the transcription starting site of AT1G79710, 684 AT4G18050, AT1G78380, AT4G24160, AT5G18480, AT4G30390, AT5G54650, AT5G16910, 685 AT5G20000, AT2G36580, AT5G56350 and AT5G17760 were amplified from Arabidopsis Col- 686 0 genomic DNA by PCR and cloned into the pGGA (plasmid Green Gate A) entry vector to 687 generate pGGA-pMarkerGene. A region of approximately -1.5 kb upstream of the transcription 688 starting site of AT1G30270 was synthetized by GENEWIZ, Inc. (South Plainfield, NJ) and 689 subsequently cloned into the pGGA entry vector as well. (Lampropoulos et al., 2013). Each 690 entry vector was then recombined with the following plasmids: pGGB-SV40-NLS, pGGC- 691 3xGFP, pGGD-RBCSt (D-F), pGGF-AlliYFP (seed coat selection cassette for transgenic seed 692 selection) and pGGZ-empty destination vector. Primers used for cloning and sequencing the 693 final constructs are listed in Table S12. All plasmids were transfected by electroporation into 694 Agrobacterium tumefaciens GV3101 strain containing the plasmid pSoup and then transformed 695 into Arabidopsis Col-0 by the floral dipping method (Clough and Bent, 1998). Transgenic 696 seeds from transformed plants were identified as those displaying a clear fluorescence signal 697 under the stereo microscope Olympus SZX18. 698 699 Pathogen inoculation and microscopy of reporter lines 700 For microscopy of the reporter line pAT5G17760:NLS-3xGFP, plants were grown as 701 previously described. Leaves of Col-0 pAT5G17760:NLS-3xGFP were infiltrated in the IN 702 area with either a mock solution (10 mM MgCl2) or different Pto strains. Pto strains expressing 703 the following effectors were used: AvrRpm1, AvrRpt2 and AvrRps4. As controls, the Pto EV 704 and Pto hrcC- strains were also used. All Pto strains were infiltrated at a wavelength of 600 705 nm (OD600) of 0.01 for microscopy imaging. For B. cinerea infection, B05.10 strain 706 was grown for 14 days in Potato Dextrose Agar (PDA) at 22°C under dark conditions. Spores 707 collected, washed in 5 mL of PDA and filtered through two layers of Miracloth (Merck 708 Millipore). Subsequently, number of spores per cm2 were counted under the microscope and 709 diluted to 1 x 105 spores per mL. For inoculation a 6 μl droplet was placed on the upper surface 710 of the 7th or 8th leaf of an adult Arabidopsis plant grown in short day conditions. A dome 711 covering the plants was placed throughout the course of B. cinerea infection. 712 713 Leaves were imaged at 16 hpi with Pto stains and 3 dpi with B. cinerea. Whole leaves were 714 photographed using a Leica DM6 microscope (Leica Microsystems) equipped with DFC365 715 FX 1.4 MP monochrome digital camera. Bright field and GFP filter pictures were taken of each 716 leaf. Confocal images were obtained using a FV1000 Olympus confocal microscope with the 717 following excitation/emission wavelengths for GFP: 488 nm/500 to 540 nm. Confocal 718 microscopy images were taken of the epidermal layer (20 Z-stacks with stack size of 1 μm) and 719 fluorescent nuclei were counted using ImageJ software. 720 721 References 722 Abu-Qare, A. W., and Abou-Donia, M. B. (2001). Biomarkers of apoptosis: Release of 723 cytochrome c, activation of caspase-3, induction of 8-hydroxy-2 ’-deoxyguanosine, 724 increased 3-nitrotyrosine, and alteration of p53 gene. Journal of Toxicology and 725 Environmental Health-Part B-Critical Reviews 4:313–332. 726 Alfano, J. R., Charkowski, A. O., Deng, W. L., Badel, J. L., Petnicki-Ocwieja, T., van 727 Dijk, K., and Collmer, A. (2000). The Pseudomonas syringae Hrp pathogenicity island 728 has a tripartite mosaic structure composed of a cluster of type III secretion genes 729 bounded by exchangeable effector and conserved effector loci that contribute to parasitic 730 fitness and pathogenicity in pl. Proc Natl Acad Sci U S A 97:4856–4861. 731 Ausubel, F. M. (2005). Are innate immune signaling pathways in plants and animals 732 conserved? Nature Immunology 6:973–979. 733 Balint-Kurti, P. (2019). The plant hypersensitive response: concepts, control and 734 consequences. Molecular Plant Pathology 20:1163–1178. 735 Barragan, A. C., Collenberg, M., Wang, J., Lee, R. R. Q., Cher, W. Y., Rabanal, F. A., 736 Ashkenazy, H., Weigel, D., and Chae, E. (2021). A Truncated Singleton NLR Causes 737 Hybrid Necrosis in Arabidopsis thaliana. Molecular Biology and Evolution 38:557–574. 738 Benjamini, Y., and Hochberg, Y. (1995). CONTROLLING THE FALSE DISCOVERY 739 RATE - A PRACTICAL AND POWERFUL APPROACH TO MULTIPLE TESTING. 740 Journal of the Royal Statistical Society Series B-Statistical Methodology 57:289–300. 741 Berardini, T. Z., Mundodi, S., Reiser, L., Huala, E., Garcia-Hernandez, M., Zhang, P. 742 F., Mueller, L. A., Yoon, J., Doyle, A., Lander, G., et al. (2004). Functional 743 annotation of the Arabidopsis genome using controlled vocabularies. Plant Physiology 744 135:745–755. 745 Berger, S., Benediktyova, Z., Matous, K., Bonfig, K., Mueller, M. J., Nedbal, L., and 746 Roitsch, T. (2007). Visualization of dynamics of plant-pathogen interaction by novel 747 combination of chlorophyll fluorescence imaging and statistical analysis: differential 748 effects of virulent and avirulent strains of P-syringae and of oxylipins on A-thaliana. 749 Journal of Experimental Botany 58:797–806. 750 Betsuyaku, S., Katou, S., Takebayashi, Y., Sakakibara, H., Nomura, N., and Fukuda, H. 751 (2018). Salicylic Acid and Jasmonic Acid Pathways are Activated in Spatially Different 752 Domains Around the Infection Site During Effector-Triggered Immunity in Arabidopsis 753 thaliana. Plant and Cell Physiology 59:8–16. 754 Bi, G., Su, M., Li, N., Liang, Y., Dang, S., Xu, J., Hu, M., Wang, J., Zou, M., Deng, Y., et 755 al. (2021). The ZAR1 resistosome is a calcium-permeable channel triggering plant 756 immune signaling. Cell 184:3528-+. 757 Bilgin, D. D., Zavala, J. A., Zhu, J., Clough, S. J., Ort, D. R., and DeLucia, E. H. (2010). 758 Biotic stress globally downregulates photosynthesis genes. Plant Cell and Environment 759 33:1597–1613. 760 Bozso, Z., Ott, P. G., Kaman-Toth, E., Bognar, G. F., Pogany, M., and Szatmari, A. 761 (2016). Overlapping Yet Response-Specific Transcriptome Alterations Characterize the 762 Nature of Tobacco-Pseudomonas syringae Interactions. Frontiers in Plant Science 7. 763 Chantarachot, T., Sorenson, R. S., Hummel, M., Ke, H., Kettenburg, A. T., Chen, D., 764 Aiyetiwa, K., Dehesh, K., Eulgem, T., Sieburth, L. E., et al. (2020). DHH1/DDX6- 765 like RNA helicases maintain ephemeral half-lives of stress-response mRNAs. Nature 766 Plants 6:675–685. 767 Chung, H. S., Koo, A. J. K., Gao, X., Jayanty, S., Thines, B., Jones, A. D., and Howe, G. 768 A. (2008). Regulation and function of Arabidopsis JASMONATE ZIM-domain genes in 769 response to wounding and herbivory. Plant Physiology 146:952–964. 770 Clough, S. J., and Bent, A. F. (1998). Floral dip: a simplified method for Agrobacterium- 771 mediated transformation of Arabidopsis thaliana. Plant Journal 16:735–743. 772 Couto, D., and Zipfel, C. (2016). Regulation of pattern recognition receptor signalling in 773 plants. Nature Reviews Immunology 16:537–552. 774 Dangl, J. L., and Jones, J. D. G. (2019). A pentangular plant inflammasome. Science (1979) 775 364:31–32. 776 Dangl, J. L., Horvath, D. M., and Staskawicz, B. J. (2013). Pivoting the Plant Immune 777 System from Dissection to Deployment. Science (1979) 341:746–751. 778 Dongus, J. A., and Parker, J. E. (2021). EDS1 signalling: At the nexus of intracellular and 779 surface receptor immunity. Current Opinion in Plant Biology 62. 780 Dorey, S., Baillieul, F., Pierrel, M. A., Saindrenan, P., Fritig, B., and Kauffmann, S. 781 (1997). Spatial and temporal induction of cell death, defense genes, and accumulation of 782 salicylic acid in tobacco leaves reacting hypersensitively to a fungal glycoprotein 783 elicitor. Molecular Plant-Microbe Interactions 10:646–655. 784 Duan, Y., Duan, S., Armstrong, M. R., Xu, J., Zheng, J., Hu, J., Chen, X., Hein, I., Li, 785 G., and Jin, L. (2020). Comparative Transcriptome Profiling Reveals Compatible and 786 Incompatible Patterns of Potato Toward Phytophthora infestans. G3-Genes Genomes 787 Genetics 10:623–634. 788 Enyedi, A. J., Yalpani, N., Silverman, P., and Raskin, I. (1992). LOCALIZATION, 789 CONJUGATION, AND FUNCTION OF SALICYLIC-ACID IN TOBACCO DURING 790 THE HYPERSENSITIVE REACTION TO TOBACCO MOSAIC-VIRUS. Proc Natl 791 Acad Sci U S A 89:2480–2484. 792 Förderer, A., Li, E., Lawson, A., Deng, Y.-N., Sun, Y., Logemann, E., Zhang, X., Wen, 793 J., Han, Z., Chang, J., et al. (2022). A wheat resistosome defines common principles of 794 immune receptor channels. bioRxiv Advance Access published 2022, 795 doi:10.1101/2022.03.23.485489. 796 Gassmann, W., Hinsch, M. E., and Staskawicz, B. J. (1999). The Arabidopsis RPS4 797 bacterial-resistance gene is a member of the TIR-NBS-LRR family of disease-resistance 798 genes. Plant Journal 20:265–277. 799 Giolai, M., Verweij, W., Lister, A., Heavens, D., Macaulay, I., and Clark, M. D. (2019). 800 Spatially resolved transcriptomics reveals plant host responses to pathogens. Plant 801 Methods 15. 802 Grant, M. R., Godiard, L., Straube, E., Ashfield, T., Lewald, J., Sattler, A., Innes, R. 803 W., and Dangl, J. L. (1995). STRUCTURE OF THE ARABIDOPSIS RPM1 GENE 804 ENABLING DUAL-SPECIFICITY DISEASE RESISTANCE. Science (1979) 269:843– 805 846. 806 Hillmer, R. A., Tsuda, K., Rallapalli, G., Asai, S., Truman, W., Papke, M. D., 807 Sakakibara, H., Jones, J. D. G., Myers, C. L., and Katagiri, F. (2017). The highly 808 buffered Arabidopsis immune signaling network conceals the functions of its 809 components. Plos Genetics 13. 810 Jacob, P., Kim, N. H., Wu, F., el Kasmr, F., Chi, Y., Walton, W. G., Furzer, O. J., 811 Lietzan, A. D., Sunil, S., Kempthorn, K., et al. (2021). Plant “helper” immune 812 receptors are Ca2+-permeable nonselective cation channels. Science (1979) 373:420-+. 813 Jones, J. D. G., and Dangl, J. L. (2006). The plant immune system. Nature 444:323–329. 814 Jones, J. D. G., Vance, R. E., and Dangl, J. L. (2016). Intracellular innate immune 815 surveillance devices in plants and animals. Science (1979) 354. 816 Kachroo, P., Burch-Smith, T. M., and Grant, M. (2021). An Emerging Role for 817 Chloroplasts in Disease and Defense. Annual Review of Phytopathology, Vol 59, 2021 818 59:423–445. 819 Kumar, L., and Futschik, M. (2007). Mfuzz: A software package for soft clustering of 820 microarray data. Bioinformation 2:5–7. 821 Lampropoulos, A., Sutikovic, Z., Wenzl, C., Maegele, I., Lohmann, J. U., and Forner, J. 822 (2013). GreenGate - A Novel, Versatile, and Efficient Cloning System for Plant 823 Transgenesis. Plos One 8. 824 Lee, M.-H., Jeon, H. S., Kim, S. H., Chung, J. H., Roppolo, D., Lee, H.-J., Cho, H. J., 825 Tobimatsu, Y., Ralph, J., and Park, O. K. (2019). Lignin-based barrier restricts 826 pathogens to the infection site and confers resistance in plants. Embo Journal 38. 827 Lewis, L. A., Polanski, K., de Torres-Zabala, M., Jayaraman, S., Bowden, L., Moore, J., 828 Penfold, C. A., Jenkins, D. J., Hill, C., Baxter, L., et al. (2015). Transcriptional 829 Dynamics Driving MAMP-Triggered Immunity and Pathogen Effector-Mediated 830 Immunosuppression in Arabidopsis Leaves Following Infection with Pseudomonas 831 syringae pv tomato DC3000. Plant Cell 27:3038–3064. 832 Littlejohn, G. R., Breen, S., Smirnoff, N., and Grant, M. (2021). Chloroplast immunity 833 illuminated. New Phytologist 229. 834 Liu, L. J., Sonbol, F. M., Huot, B., Gu, Y. N., Withers, J., Mwimba, M., Yao, J., He, S. 835 Y., and Dong, X. N. (2016). Salicylic acid receptors activate jasmonic acid signalling 836 through a non-canonical pathway to promote effector-triggered immunity. Nature 837 Communications 7:10. 838 Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and 839 dispersion for RNA-seq data with DESeq2. Genome Biology 15. 840 Love, M. I., Anders, S., Kim, V., and Huber, W. (2015). RNA-Seq workflow: gene-level 841 exploratory analysis and differential expression. F1000Res 4:1070. 842 Lu, Y., and Yao, J. (2018). Chloroplasts at the Crossroad of Photosynthesis, Pathogen 843 Infection and Plant Defense. International Journal of Molecular Sciences 19. 844 Lu, J. v, Chen, H. C., and Walsh, C. M. (2014). Necroptotic signaling in adaptive and 845 innate immunity. Seminars in Cell & Developmental Biology 35:33–39. 846 Lukan, T., Pompe-Novak, M., Baebler, S., Tusek-Znidaric, M., Kladnik, A., Kriznik, 847 M., Blejec, A., Zagorscak, M., Stare, K., Dusak, B., et al. (2020). Precision 848 transcriptomics of viral foci reveals the spatial regulation of immune-signaling genes 849 and identifiesRBOHDas an important player in the incompatible interaction between 850 potato virus Y and potato. Plant Journal 104:645–661. 851 Ma, S., Lapin, D., Liu, L., Sun, Y., Song, W., Zhang, X., Logemann, E., Yu, D., Wang, 852 J., Jirschitzka, J., et al. (2020). Direct pathogen-induced assembly of an NLR immune 853 receptor complex to form a holoenzyme. Science (1979) 370:1184-+. 854 Mackey, D., Holt, B. F., Wiig, A., and Dangl, J. L. (2002). RIN4 interacts with 855 Pseudomonas syringae type III effector molecules and is required for RPM1-mediated 856 resistance in Arabidopsis. Cell 108:743–754. 857 Mackey, D., Belkhadir, Y., Alonso, J. M., Ecker, J. R., and Dangl, J. L. (2003). 858 Arabidopsis RIN4 is a target of the type III virulence effector AvrRpt2 and modulates 859 RPS2-mediated resistance. Cell 112:379–389. 860 Martin, R., Qi, T., Zhang, H., Liu, F., King, M., Toth, C., Nogales, E., and Staskawicz, 861 B. J. (2020). Structure of the activated ROQ1 resistosome directly recognizing the 862 pathogen effector XopQ. Science (1979) 370:1185-+. 863 Mine, A., Seyfferth, C., Kracher, B., Berens, M. L., Becker, D., and Tsuda, K. (2018). 864 The Defense Phytohormone Signaling Network Enables Rapid, High-Amplitude 865 Transcriptional Reprogramming during Effector-Triggered Immunity. Plant Cell 866 30:1199–1219. 867 Misas-Villamil, J. C., Kolodziejek, I., Crabill, E., Kaschani, F., Niessen, S., Shindo, T., 868 Kaiser, M., Alfano, J. R., and van der Hoorn, R. A. L. (2013). Pseudomonas syringae 869 pv. syringae Uses Proteasome Inhibitor Syringolin A to Colonize from Wound Infection 870 Sites. Plos Pathogens 9. 871 Muckenschnabel, I., Goodman, B. A., Williamson, B., Lyon, G. D., and Deighton, N. 872 (2002). Infection of leaves of Arabidopsis thaliana by Botrytis cinerea: Changes in 873 ascorbic acid, free radicals and lipid peroxidation products. Journal of Experimental 874 Botany 53:207–214. 875 Nagata, S., and Tanaka, M. (2017). Programmed cell death and the immune system. Nature 876 Reviews Immunology 17:333–340. 877 Narusaka, M., Shirasu, K., Noutoshi, Y., Kubo, Y., Shiraishi, T., Iwabuchi, M., and 878 Narusaka, Y. (2009). RRS1 and RPS4 provide a dual Resistance-gene system against 879 fungal and bacterial pathogens. Plant Journal 60:218–226. 880 Navarro, L., Zipfel, C., Rowland, O., Keller, I., Robatzek, S., Boller, T., and Jones, J. D. 881 G. (2004). The transcriptional innate immune response to flg22. interplay and overlap 882 with Avr gene-dependent defense responses and bacterial pathogenesis. Plant 883 Physiology 135:1113–1128. 884 Ngou, B. P. M., Jones, J. D. G., and Ding, P. (2021a). Plant immune networks. Trends 885 Plant Sci Advance Access published 2021, doi:10.1016/j.tplants.2021.08.012. 886 Ngou, B. P. M., Ahn, H. K., Ding, P., and Jones, J. D. G. (2021b). Mutual potentiation of 887 plant immunity by cell-surface and intracellular receptors. Nature 2021 592:7852 888 592:110–115. 889 Olsen, J. v, Blagoev, B., Gnad, F., Macek, B., Kumar, C., Mortensen, P., and Mann, M. 890 (2006). Global, in vivo, and site-specific phosphorylation dynamics in signaling 891 networks. Cell 127:635–648. 892 Olvera-Carrillo, Y., van Bel, M., van Hautegem, T., Fendrych, M., Huysmans, M., 893 Simaskova, M., van Durme, M., Buscaill, P., Rivas, S., Coll, N. S., et al. (2015). A 894 Conserved Core of Programmed Cell Death Indicator Genes Discriminates 895 Developmentally and Environmentally Induced Programmed Cell Death in Plants. Plant 896 Physiology 169:2684–2699. 897 Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., and Kingsford, C. (2017). Salmon 898 provides fast and bias-aware quantification of transcript expression. Nature Methods 899 14:417-+. 900 Pitsili, E., Phukan, U. J., and Coll, N. S. (2020). Cell Death in Plant Immunity. Cold Spring 901 Harbor Perspectives in Biology 12. 902 Pruitt, R. N., Locci, F., Wanke, F., Zhang, L., Saile, S. C., Joe, A., Karelina, D., Hua, C., 903 Frohlich, K., Wan, W.-L., et al. (2021). The EDS1-PAD4-ADR1 node mediates 904 Arabidopsis pattern-triggered immunity. Nature Advance Access published 2021, 905 doi:10.1038/s41586-021-03829-0. 906 RStudio | Open source & professional software for data science teams (2021).Advance 907 Access published 2021. 908 Salguero-Linares, J., and Coll, N. S. (2019). Plant proteases in the control of the 909 hypersensitive response. Journal of Experimental Botany 70:2087–2095. 910 Savage, Z., Duggan, C., Toufexi, A., Pandey, P., Liang, Y., Eugenia Segretin, M., Yuen, 911 L. H., Gaboriau, D. C. A., Leary, A. Y., Tumtas, Y., et al. (2021). Chloroplasts alter 912 their morphology and accumulate at the pathogen interface during infection by 913 Phytophthora infestans. Plant Journal Advance Access published 2021, 914 doi:10.1111/tpj.15416. 915 Savatin, D. v, Gramegna, G., Modesti, V., and Cervone, F. (2014). Wounding in the plant 916 tissue: the defense of a dangerous passage. Frontiers in Plant Science 5. 917 Schreiber, U. (2004). Pulse-amplitude-modulation (PAM) fluorometry and saturation pulse 918 method: An overview. Chlorophyll a Fluoerescence: Signature of Photosynthesis 919 19:279-+. 920 Schreiber, U., Klughammer, C., and Kolbowski, J. (2012). Assessment of wavelength- 921 dependent parameters of photosynthetic electron transport with a new type of multi- 922 color PAM chlorophyll fluorometer. Photosynthesis Research 113:127–144. 923 Soneson, C., Love, M. I., and Robinson, M. D. (2015). Differential analyses for RNA-seq: 924 transcript-level estimates improve gene-level inferences. F1000Res 4:1521. 925 Spoel, S. H., Johnson, J. S., and Dong, X. (2007). Regulation of tradeoffs between plant 926 defenses against pathogens with different lifestyles. Proc Natl Acad Sci U S A 927 104:18842–18847. 928 Tao, Y., Xie, Z. Y., Chen, W. Q., Glazebrook, J., Chang, H. S., Han, B., Zhu, T., Zou, G. 929 Z., and Katagiri, F. (2003). Quantitative nature of Arabidopsis responses during 930 compatible and incompatible interactions with the bacterial pathogen Pseudomonas 931 syringae. Plant Cell 15:317–330. 932 Tian, H., Wu, Z., Chen, S., Ao, K., Huang, W., Yaghmaiean, H., Sun, T., Xu, F., Zhang, 933 Y., Wang, S., et al. (2021). Activation of TIR signalling boosts pattern-triggered 934 immunity. Nature Advance Access published 2021, doi:10.1038/s41586-021-03987-1. 935 Ustun, S., Sheikh, A., Gimenez-Ibanez, S., Jones, A., Ntoukakis, V., and Bornke, F. 936 (2016). The Proteasome Acts as a Hub for Plant Immunity and Is Targeted by 937 Pseudomonas Type III Effectors. Plant Physiology 172:1941–1958. 938 Ustun, S., Hafren, A., Liu, Q. S., Marshall, R. S., Minina, E. A., Bozhkov, P. v, Vierstra, 939 R. D., and Hofius, D. (2018). Bacteria Exploit Autophagy for Proteasome Degradation 940 and Enhanced Virulence in Plants. Plant Cell 30:668–685. 941 van Doorn, W. G. (2011). Classes of programmed cell death in plants, compared to those in 942 animals. Journal of Experimental Botany 62:4749–4761. 943 van Doorn, W. G., Beers, E. P., Dangl, J. L., Franklin-Tong, V. E., Gallois, P., Hara- 944 Nishimura, I., Jones, A. M., Kawai-Yamada, M., Lam, E., Mundy, J., et al. (2011). 945 Morphological classification of plant cell deaths. Cell Death and Differentiation 946 18:1241–1246. 947 Vega-Munoz, I., Duran-Flores, D., Fernandez-Fernandez, A. D., Heyman, J., Ritter, A., 948 and Stael, S. (2020). Breaking Bad News: Dynamic Molecular Mechanisms of Wound 949 Response in Plants. Frontiers in Plant Science 11. 950 Wang, J., Hu, M., Wang, J., Qi, J., Han, Z., Wang, G., Qi, Y., Wang, H.-W., Zhou, J.- 951 M., and Chai, J. (2019a). Reconstitution and structure of a plant NLR resistosome 952 conferring immunity. Science (1979) 364:44-+. 953 Wang, J. Z., Wang, J., Hu, M. J., Wu, S., Qi, J. F., Wang, G. X., Han, Z. F., Qi, Y. J., 954 Gao, N., Wang, H. W., et al. (2019b). Ligand-triggered allosteric ADP release primes a 955 plant NLR complex. Science (1979) 364:43-+. 956 Ward, T., Cummings, J., Dean, E., Greystoke, A., Hou, J. M., Backen, A., Ranson, M., 957 and Dive, C. (2008). Biomarkers of apoptosis. British Journal of Cancer 99:841–846. 958 Wu, C. H., Abd-El-Haliem, A., Bozkurt, T. O., Belhaj, K., Terauchi, R., Vossen, J. H., 959 and Kamoun, S. (2017). NLR network mediates immunity to diverse plant pathogens. 960 Proc Natl Acad Sci U S A 114:8113–8118. 961 Wu, C.-H., Derevnina, L., and Kamoun, S. (2018). Receptor networks underpin plant 962 immunity. Science (1979) 360:1300–1301. 963 Yang, L., Chen, X., Wang, Z., Sun, Q., Hong, A., Zhang, A., Zhong, X., and Hua, J. 964 (2020). HOS15 and HDA9 negatively regulate immunity through histone deacetylation 965 of intracellular immune receptor NLR genes in Arabidopsis. New Phytologist 226:507– 966 522. 967 Yu, D., Song, W., Yong, E., Tan, J., Liu, L., Cao, Y., Jirschitzka, J., Li, E., Logemann, 968 E., Xu, C., et al. (2021). TIR domains of plant immune receptors are 2′,3′-cAMP/cGMP 969 synthetases mediating cell death. bioRxiv Advance Access published November 10, 970 2021, doi:10.1101/2021.11.09.467869. 971 Yuan, M., Jiang, Z., Bi, G., Nomura, K., Liu, M., Wang, Y., Cai, B., Zhou, J.-M., Yang 972 He, S., Xin, X.-F., et al. (2021a). Pattern-recognition receptors are required for NLR- 973 mediated plant immunity Check for updates. Nature 592:105. 974 Yuan, M., Ngou, B. P. M., Ding, P., and Xin, X.-F. (2021b). PTI-ETI crosstalk: an 975 integrative view of plant immunity. Curr Opin Plant Biol 62:102030. 976 Zheng, X., Zhou, M., Yoo, H., Pruneda-Paz, J. L., Spivey, N. W., Kay, S. A., and Dong, 977 X. (2015). Spatial and temporal regulation of biosynthesis of the plant immune signal 978 salicylic acid. Proc Natl Acad Sci U S A 112:9166–9173. 979 980 Author Contribution Statement 981 JS-L designed and performed experiments, analyzed and interpreted data and wrote the 982 manuscript 983 IS designed and performed experiments and analyzed and interpreted data and helped writing 984 the manuscript 985 NR-S designed and performed experiments, analyzed and interpreted data and helped writing 986 the manuscript 987 MS performed experiments 988 UP performed experiments 989 VMG performed analysis and interpreted data 990 MB-F performed analysis and interpreted data 991 MV interpreted data and helped writing the manuscript 992 DR performed experiments, analyzed and interpreted data, and helped writing the manuscript. 993 NSC conceptualized the research, designed the experiments, interpreted data and wrote the 994 manuscript. 995 996 Acknowledgements 997 The authors would like to thank Susana Rivas, who conceived and initiated the project, but 998 declined to be author on the manuscript. Likewise, we thank Susana Rivas’s team for their help 999 with the preliminary experiments and the plant tissue harvest for the RNA-Seq. We thank 1000 Sebastien Carrère from the Bioinformatics facility at the LIPM, for his bioinformatics 1001 preliminary analysis. We also thank Simon Stael (VIB) for helpful comments and inspiring 1002 discussions and all members from the Bacterial plant diseases and cell death lab for their 1003 insights and suggestions. We thank José Luis Riechman (CRAG) and Miguel Ángel Moreno- 1004 Risueño for help with the analysis, Antoni Garcia-Molina for help with B. cinerea infections 1005 and Montse Amenós for help with microscopy. We would like to thank Kenichi Tsuda for 1006 sharing his RNA-seq data of previously published transcriptomic studies (Mine et al., 2018) 1007 and Ignacio Rubio-Somoza for providing us with the green gate plasmid pGGD-RBCSt (D-F). 1008 1009 Conflict of Interest Statement 1010 The authors declare no conflict of interest. 1011 1012 Ethics Statement 1013 The present study did not require ethical approval. 1014 1015 Funding Statement 1016 Research at CRAG was funded with grants PID2019-108595RB-I00 funded by MCIN/AEI/ 1017 10.13039/501100011033 and AGL2016-78002-R funded by MCIN/AEI/ 1018 10.13039/501100011033 and by “ERDF A way of making Europe” (NSC, MV), fellowship 1019 PID2019-108595RB-I00 funded by Spanish MCIN/AEI/ 10.13039/501100011033 (NSC, 1020 MV) and fellowships BES-2017-080210 funded by MCIN/AEI/ 10.13039/501100011033 and 1021 by “ESF Investing in your future” (JS-L) and FPU19/03778 funded by MU (o Ministerio de 1022 Universidades) (NR-S); and through the “Severo Ochoa Programme for Centres of Excellence 1023 in R&D” (SEV-2015-0533 and CEX2019-000902-S funded by MCIN/AEI/ 1024 10.13039/501100011033) and by the CERCA Programme / Generalitat de Catalunya. Work at 1025 the LIPM was supported by the INRA SPE department (AAP2014), the Région Midi-Pyrénees 1026 (grant 13050322) and the French Laboratory of Excellence project “TULIP” (ANR-10-LABX- 1027 41; ANR-11-IDEX-0002-02). IS was supported by an AgreenSkills fellowship within the EU 1028 Marie-Curie FP7 COFUND People Programme (grant agreement no. 267196). 1029 1030 Data Availability Statement 1031 RNA-seq raw and processed data generated in this study can be found in GEO (GSE198022). 1032 All code used for analysis can be found at 1033 https://doi.org/10.34810/data174 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 Figure legends Figure 1. HR in plants can be spatio-temporally dissected. (A) Experimental design of the study. A limited area (3-4 mm) at the side edge of four-week-old Arabidopsis thaliana Col-0 leaves was syringe-infiltrated with either Pto AvrRpm1 at 2.5*107 cfu/ml (INFECTED) or a 10 mM MgCl2 solution (MOCK) and samples were collected at 5 different time points after infection: 0, 1, 2, 4 and 6 hpi. Upon infiltration, the edge of the infiltrated area was marked, and the total area infiltrated designated as “IN”. A 1 mm buffer zone right next to the IN zone ensured proper separation between the IN and “OUT” area, which was the parallel region that expanded from the edge of the buffer zone to 1-2 mm towards the vein. Three biological replicates per area, treatment and time point were collected and subjected for RNA-seq analysis. (B) Analysis of macroscopic cell death upon infection with either Pto AvrRpm1 or 10 mM MgCl2 solution. Leaves were infected as described in (a) and subsequently stained with trypan blue. Scale bar 3 mm (C) Representative images of mock or Pto AvrRm1-treated plants subjected to pulse-amplitude modulated (PAM) chlorophyll fluorescence measurement to monitor photosynthesis. Scale bar 3 mm. Photosynthetic efficiency (Fv/Fm ratio) and electron transport rate (ETR) were measured in the infiltrated area (IN) and the neighboring tissue (OUT). Measurements were taken at 0, 1, 2, 4 and 6 hpi. Results are representative of 6 different measurements of each tissue area from 6 different plants. Letters indicate statistically significant differences in either Fv/Fm ratio or ETR values following a two-way ANOVA with Tukey´s HSD test (α = 0.05). Exact p values are provided in Table S5. Figure 2. Spatio-temporal dynamics of the transcriptome reveal time and zone-dependent gene expression signatures upon infection. (A) Differentially expressed genes (FDR < 0.05 and |log2FC| > 2) in Pto AvrRpm1-infected plants compared to mock-treated plants at each time point at the IN (left) and OUT (right) areas. Red color denotes upregulated genes whereas blue color indicates downregulated genes. Yellow color shows genes with FDR < 0,05 but |log2FC| < 2, whereas grey indicates genes not complying with neither FDR nor log2FC criteria. (B-C) Genes exclusively upregulated (FDR < 0.05 and log2FC>2) at either IN or OUT areas of infection at 4 and 6 hpi. (B) Venn diagrams shows sizes of gene sets that are upregulated (FDR < 0.05 and log2FC > 2) upon bacterial infection at 4 and/or 6 hpi at either IN, OUT or both areas. (C) GO terms representing enriched biological processes derived from genes exclusively upregulated at either IN or OUT areas at 4 and/or 6 hpi. The most specific term from each family term provided by PANTHER was plotted along with their corresponding gene number, fold enrichment (FE) and FDR (Bonferroni Correction for multiple testing) represented as log10. Only GO Terms with a FE above 2 and FDR below 0.05 were plotted. Figure 3. Gene expression profile clustering reveals three distinctive expression patterns at the IN and OUT areas of infection. Non-overlapping clusters derived from Pto AvrRpm1- and mock-treated plants for IN (A) and OUT (B) areas. Standardized expression to Z-scores (Y-axis) is calculated by subtracting the mean and normalizing to standard deviation. The trajectory that defines the overall expression profile of each cluster through the course of the infection is shown in red for Pto AvrRpm1-treated plants. Genes derived from Pto AvrRpm1- treated samples were re-clustered for mock-treated samples and their trajectories are represented in grey. Since the expression profile of these genes in mock-treated samples was very distinct among the overall number of genes, they were divided into two sub-clusters represented either in dotted or dashed grey lines. The number of genes that constitute each cluster is indicated below each cluster. Genes comprising each cluster along with their MSV can be found in Table S6-S7. Figure 4. Identification of HR markers specific for the IN area of infection. (A) Schematic representation of the sequence of filters applied to identify indicators. Four filters were concatenated considering the three variables of our experimental design: time, treatment and tissue area. Briefly, in the first filter, we selected genes differentially upregulated from 0 to 4/6 hpi (FDR < 0.05 and log2FC > 2) at the IN area (colored in red) upon bacterial infection. From the genes that passed this first filter, we selected those that were exclusively upregulated (FDR < 0.05 and log2FC > 2) due to bacterial infection at the IN area at 4/6 hpi. Subsequently, from the genes that made it into the third filter, we selected those that were not highly upregulated in the OUT area (colored in blue) upon bacterial infection at 4/6 hpi (FDR < 0.05 and log2FC < 1). Finally, we applied a fourth filter to discard genes that could potentially be basally upregulated at the OUT area upon pathogen treatment at 4/6 hpi (FDR < 0.05 and log2FC > 2). The starting number of genes and the genes passing the different filtering criteria are indicated. (B) RT-qPCR and RNA-seq expression profiles of marker genes that behave as bona fide HR indicators. Relative expression levels to the housekeeping gene EIF4a were represented as fold enrichment between 4/6 and 0 hpi. Error bars represent standard error of the mean from three independent experiments. Letters indicate statistically significant differences between treatments following one-way ANOVA with Tukey´s HSD test (α = 0.05) performed independently at IN and OUT. NS (non-significant after one-way ANOVA). Exact p values are provided in Table S5. (C) List of HR indicators along with their gene ID, gene name and description. igure 5. AT5G17760 encodes an AAA-ATPase and is a reliable HR indicator specifically induced at the IN area by activation of different classes of NLR receptors. (A) Representative images of trypan blue-stained leaves, epifluorescence microscope and confocal microscope from pAT5G17760:NLS-3xGFP Arabidopsis transgenics. A small region of 4- week-old pAT5G17760::NLS-3xGFP leaves was syringe-infiltrated with Pto expressing the effectors AvrRpm1, AvrRpt2 or AvrRps4 at 1*107 colony-forming units (CFU)/ml (O.D600 = 0.01). Besides mock treatment, the non-cell death-causing bacterial strains Pto DC3000 EV and Pto DC3000 hrcC- were included as negative controls. Images were taken 16 hpi. Scale bar 3 mm. Images were taken 16 hpi on a Leica DM6 microscope and a confocal microscope prior to trypan blue staining. Scale bar 3 mm. Expression of pAT5G17760 is detected as green dots corresponding to nuclei with positive GFP signal. Scale bar 100 µm. A representative close-up magnified image of a Pto AvrRpm1-infected leaf expressing pAT5G17760:NLS- 3xGFP at 16 hpi is shown. Scale bar 3 mm. (B) Quantification of fluorescent nuclei from confocal microscopy pictures in (A). Nuclei count was performed using ImageJ software. Data is representative of three independent experiments each one of them containing 4 leaves. Letters indicate statistically significant differences in number of nuclei following one-way ANOVA with Tukey´s HSD test (α = 0.05). Exact p values are provided in Table S5. (C) Activation of pAT5G17760 occurred upon drop inoculation infection with the necrotrophic pathogen B. cinerea or a mock solution. Four- to 5-week-old leaves from pAT5G17760::3xGFP transgenics were drop inoculated with B. cinerea at concentration of 1x105 spores per mL. Images in left panels represent trypan blue stained leaves and right panels represent leaves imaged under the epifluoresecent microscope at 3 dpi. Scale bars 3 mm. Figure 1 Click here to access/download;Figure;Figure 1.TIF Figure 2 Click here to access/download;Figure;Figure 2.TIF Figure 3 Click here to access/download;Figure;Figure 3.TIF Figure 4 Click here to access/download;Figure;Figure 4.TIF Figure 5 Click here to access/download;Figure;Figure 5.TIF Supplemental Information Robust transcriptional indicators of immune cell death revealed by spatio- temporal transcriptome analyses Jose Salguero-Linaresa,#, Irene Serranob, ,#,†, Nerea Ruiz-Solania, Marta Salas-Gómeza, Ujjal Jyoti Phukana, Victor Manuel Gonzáleza, Martí Bernardo-Fauraa, Marc Vallsa,b, David Rengelb,c,¥,§,*, Nuria S. Colla,d,§,* Supplemental Information Supplementary figures and supplementary tables Figure S1. Principal component analysis (PCA) from the RNA seq-data. Circles represent mock-treated plants and triangles represent Pto AvrRpm1-infected plants. Different colors are assigned for each time point. (A) PCA comprising all data sets in our study (IN and OUT samples together). (B) PCA with IN and OUT data sets separated in order to ease visualization of the data. Figure S2. GO term enrichment analysis of upregulated and downregulated genes at each time after infection at the IN (A-B) and OUT (C) areas. The most specific term from each family term provided by PANTHER was plotted along with their corresponding gene number, fold enrichment and adj p value (Bonferroni Correction for multiple testing) represented as log10. Only GO terms with a fold enrichment above 2 and adj p value below 0.05 were plotted. Figure S3. The majority of differentially expressed genes at both IN and OUT are specific to 4 and 6 hpi. Venn diagrams showing sizes of gene sets that are differentially expressed (red: upregulated and blue: downregulated) at IN (A) or OUT (B) at each time point. Figure S4. Heatmap representing differential expression of genes exclusively upregulated at 4 and/or 6 hpi at the OUT area (log2FC > 2 and BTH <0.05) throughout the course of the infection (0,1,2,4 and 6 hpi) at IN and OUT areas. Figure S5. RNA-seq expression profiles of JA responsive genes exclusively upregulated at the OUT area upon Pto AvrRpm1 infection. Gene expression of genes from Pto-AvrRpm1 or mock- infected plants is represented as DESeq2 pseudocounts. JAL35, Jacalin-related lectin 35; CYT1, Mannose-1-phosphate guanylyltransferase 1; 4CLL5, 4-coumarate--CoA ligase-like 5; TIFY7, Protein TIFY 7; CYP74A, Allene oxide synthase, chloroplastic; RGL3,, DELLA protein RGL3; TIFY6B, Protein TIFY 6B; TIFY10B, Protein TIFY 10B; JAR1, Jasmonoyl--L-amino acid synthetase JAR1; NPF6.2, Protein NRT1/ PTR FAMILY 6.2 Figure S6. GO terms representing enriched biological processes derived from each cluster in Pto AvrRpm1-treated plants. GO term enrichment analysis was performed on those genes that had a membership score value (MSV) above or equal to 0.7 (see Materials and Methods). The most specific term from each family provided by PANTHER was plotted along with their corresponding gene number, fold enrichment (FE) and adj p value (Bonferroni Correction for multiple testing) represented as log10. Only GO Terms with a FE above 2 and adj p value below 0.05 were plotted. Enriched GO terms from cluster I (2,937 genes; MSV > 0.7  1069 genes), cluster II (4,183 genes; MSV > 0.7  2613 genes) and cluster III (6,428 genes; MSV > 0.7  4885 genes) at the IN area (A) in Pto AvrRpm1-treated plants were predominantly linked to processes related to immunity, protein turnover and photosynthesis, respectively. At the OUT area (B), enriched GO terms from cluster I (1,552 genes; MS > 0.7  747 genes) and II (1,100 genes; MS > 0.7  184) suggest the importance of processes related to hormonal regulation in by-stander cells, whereas genes comprising cluster III (925 genes; MS > 0.7  181 genes) infer that photosynthesis and rearrangements in the chloroplast occur similarly compared to mock-treated samples at the OUT area Figure S7. GO terms representing enriched biological processes derived from each sub-cluster in mock-treated plants at the IN and OUT areas. From each cluster belonging to mock-treated samples, GO term enrichment analysis was performed on those genes that had a membership score value (MSV) above or equal to 0.7 at the IN (a) and OUT areas (b). The most specific term from each family term provided by PANTHER was plotted along with their corresponding gene number, fold enrichment and adj p value (Bonferroni Correction for multiple testing) represented as log10. Only GO Terms with a fold enrichment above 2 and adj p value below 0.05 were plotted. (A) Sub-cluster 1.1 (638 genes; MSV >= 0.7  467 genes), sub-cluster 1.2 (2299 genes; MSV >= 0.7  1942 genes), sub-cluster 2.1 (2570 genes; MSV >= 0.7  1573 genes), sub-cluster 2.2 (1613 genes; MSV >= 0.7  649 genes), sub-cluster 3.1 (3172 genes; MSV >= 0.7  2391 genes), sub-cluster 3.2 (3256 genes; MSV >= 0.7  2557 genes). (B) Sub-cluster 1.1 (850 genes; MSV >= 0.7  319 genes), sub-cluster 1.2 (702 genes; MSV >= 0.7  183 genes), sub-cluster 2.1 (453 genes; MSV >= 0.7  286 genes), sub-cluster 2.2 (647 genes; MSV >= 0.7  389 genes), sub-cluster 3.1 (612 genes; MSV >= 0.7  555 genes), sub-cluster 3.2 (313 genes; MSV >= 0.7  257 genes). Figure S8. Clusters 4 (1,174 genes; MSV >= 0.7  57 genes ) and 5 (961 genes; MSV >= 0.7  314 genes) from Pto AvrRpm1-treated plants at the OUT area share similar expression profiles and do not contain any relevant enriched GO terms associated with biological processes, possibly due to low gene number. Figure S9. List of in silico HR indicators obtained after filtering at 4 and 6 hpi. (A) Briefly, we firstly selected genes that were upregulated (log2FC > 2) after Pto AvrRpm1 infection at 4 or 6 hpi vs 0 hpi. From the genes that complied with this first filter, we selected those that were specifically upregulated in Pto AvrRpm1-infected vs mock-inoculated samples at 4 or 6 hpi (log2FC >2). From the genes that complied these criteria, we kept those with a log2FC <1 at the OUT area in Pto AvrRpm1-infected vs mock-inoculated samples at 4 or 6 hpi. Finally, from the genes that met those three criteria, we kept those that were differentially upregulated at the IN area compared to the OUT area in Pto AvrRpm1-infected plants. (B) Log2FCs resulting from pairwise comparisons in the 1st, 2nd, 3rd and 4th filters applied are indicated for each gene marker along with its corresponding gene description. Figure S10. RNA-seq expression profiles of 4 (A) and 6 (B) hour candidate HR indicators at the IN and OUT areas of infection. Gene expression of genes from Pto-AvrRpm1 or mock- infected plants is represented as DESeq2 pseudocounts. Figure S11. RT-qPCR of 4- and 6-hour transcriptional HR indicators at IN and OUT areas upon treatment with either mock, Pto AvrRpm1 or Pto DC3000 EV. Relative expression levels to the housekeeping gene EIF4a were represented as fold induction between 4 (A) or 6 (B) and 0 hpi. Error bars represent standard error of the mean from three independent experiments. Letters indicate statistically significant differences between treatments following one-way ANOVA with Tukey´s HSD test (α = 0.05) performed independently at IN and OUT. NS (non- significant after one-way ANOVA). Exact p values are provided in Table S5. Figure S12. Time course imaging of pAT5G17760:NLS-3xGFP Arabidopsis transgenic leaves infected with Pto AvrRpm1 (A), Pto DC3000 EV (B) or mock solution (10 mM MgCl2) (C). A small region of 4-week-old pAT5G17760::NLS-3xGFP leaves was syringe-infiltrated with Pto strains at 1*107 colony-forming units (CFU)/ml (O.D600 = 0.01). Fluorescent microscopy images were taken at 0, 2, 4, 8, 12, 16 and 24 hpi (right panels). Afterwards, leaves were subjected to trypan blue staining (left panels). Scale bar 3 mm. Figure S13. Activation of pAT5G17760 in the syringe-infiltrated area occurred in several independent pAT5G17760::3xGFP transgenic lines. Leaves of Arabidopsis transgenics in the T2 generation were syringe infiltrated with Pto AvrRpm1 (A), Pto DC3000 EV (B) at 1*107 colony-forming units (CFU)/ml (O.D600 = 0.01) and imaged at 16 hpi. Mock solution was used as a control (C). Images in left panels are leaves stained with trypan blue whereas images in right panels are leaves under the epifluorescence microscope. Scale bars 3 mm. Figure S14. Disease resistance and cell death triggered by avirulent Pto AvrRpm1 strain is not compromised in Arabidopsis mutant lacking AT5G17760. (A) Scheme of AT517760 gene indicating the position of the T-DNA insertion in GK-59F04 mutant line (left panel) and RT-qPCR of two regions (F1-R1 and F2-R2) of exon 2 in Col-0 and GK-59F04 plants. RT- qPCR data is represented as relative expression levels of AT5G17760 to the housekeeping gene EIF4a (right panel). Error bars represent standard error of the mean from four biological replicates. Letters indicate statistically significant differences between treatments following a Welch Two Sample t-test. Exact p values are provided in Table S5. (B-C) Four to 5 week-old Col-0, at5g17760 and rpm1-3 plants were syringe-infiltrated with Pto DC3000 AvrRpm1 at O.D600=0.05 for electrolyte leakage (b) and O.D600=0.001 for bacterial growth assays (C), rpm1-3 mutant is used as a negative control since it is defective in the cognate NLR that recognizes the effector AvrRpm1. (B) Conductivity measurements of electrolyte leakage from dying cells were recorded at 0, 4, 6, 8 and 10 hpi. Dots represent data from 3 biological replicates (represented in different colors) consisting of 4 technical replicates each with 2 leaf discs measured per replicate. Letters indicate statistically significant differences between genotypes following one-way ANOVA with Tukey´s HSD test performed at each time point. Exact p values are provided in Table S5. (C) Bacterial growth at 0 and 3 days post-infection (dpi) was measured in Col-0, at5g17760 and rpm1-3. Dots represent bacterial CFU (colony- forming units) per cm2 from 2 biological replicates (represented in different colors) consisting of 4 technical replicates each with 2 leaf discs measured per replicate. Letters indicate statistically significant differences between genotypes following one-way ANOVA with Tukey´s HSD test performed at 0 and 3 days post infection. Exact p values are provided in Table S5. Figure S15. RNA-seq expression profiles of genes involved in lignin biosynthesis. (A) Gene expression of genes from Pto-AvrRpm1 or mock-infected plants is represented as DESeq2 pseudocounts. (B) Scheme of lignin biosynthesis in plants. Black arrow indicates the canonical lignin biosynthesis in plants. Bold font indicates enzymes involved in the different steps of the pathway. PAL, phenylalanine ammonia-lyase; C4H, cinnamate 4-hydroxylase; 4CL, 4-coumarate: CoA ligase; HCT, quinateshikimate p-hydroxycinnamoyltransferase; C3′H, p-coumaroylshikimate 3′-hydroxylase; CCoAOMT, caffeoyl-CoA O- methyltransferase; CCR, cinnamoyl-CoAreductase; F5H, ferulate 5-hydroxylase; CAD, cinnamyl alcohol dehydrogenase; COMT, caffeic acid O-methyltransferase; CSE, caffeoyl shikimate esterase; PRX, peroxidase; LAC, laccase (Adapted from Meng Chie et al., 2018) Figure S16. Transcriptional regulation of HR markers found in this study compared to RNA-seq data sets from plants undergoing ETI and autoimmunity. (A) HR markers found in this study were searched in Mine et al., 2018. In pairwise comparisons between infection with ETI-causing bacteria (Pto AvrRpm1 and Pto AvrRpt2) and mock, only genes with high statistical confidence (q value ≤ 0,01) in at least one time point were plotted on a heatmap indicating Log2FC for the different times tested in their study. (B) HR markers found in this study were searched in RNA-seq data sets of Arabidopsis hos15-4 (Yang et al.,2019), rh6812 autoimmune plants (Chantarachot et al., 2020) and Cdm-0 x TueScha-9 F1 hybrids (Barragan et al., 2020). Genes with high statistical confidance in their data sets (FDR <0.05) were plotted on a heatmap indicating Log2FC between expression of WT (Col-0) and autoimmune plants. Supplementary tables/dataset legends. Table S1 (Associated to Figure 2A). List of differentially expressed genes upon Pto AvrRpm1 infection at each time point and tissue area. Table S2 (Associated to Figure S2). List of genes constituting each GO term in Figure S2. GO term enrichment analysis of upregulated and downregulated genes at either IN our OUT areas. Only those GO terms exhibiting an FDR < 0.05 after Bonferroni Correction for multiple testing and a fold enrichment above 2 are shown. Table S3 (Associated to Figure 2B). List of genes that are upregulated upon Pto AvrRpm1 at 4 and 6 hpi exclusively at IN, both at IN and OUT or exclusively at OUT. Table S4 (Associated to Figure 2C). List of genes constituting each GO term in Figure 2C.GO term enrichment analysis of genes that are exclusively upregulated at either the IN or OUT area upon Pto AvrRpm1 infection. Only those GO terms exhibiting an FDR < 0.05 after Bonferroni Correction for multiple testing and a fold enrichment above 2 are shown. Table S5 (Associated to Figure 1C, Figure 4, Figure 5, Figure S11 and Figure S14). Tukey HSD p-values and Welch two sample t-test p-values obtained from statistical tests applied in the study. Table S6 (Associated to Figure 3A). List of genes comprising each cluster derived from Pto AvrRpm1 and mock-treated plants along with their corresponding MSV at IN. Table S7 (Associated to Figure 3B). List of genes comprising each cluster derived from Pto AvrRpm1 and mock-treated plants along with their corresponding MSV at OUT. Table S8 (Associated to Figure S6A). List of genes constituting each GO term in Figure S6A. GO term enrichment analysis of genes from clusters of Pto AvrRpm1-inoculated plants at the IN area with a MSV of 0.7 or above. Only those GO terms exhibiting an FDR < 0.05 after Bonferroni Correction for multiple testing and a fold enrichment above 2 are shown. Table S9 (Associated to Figure S6B). List of genes constituting each GO term in Figure S6B. GO term enrichment analysis of genes from clusters of Pto AvrRpm1-inoculated plants at the OUT area with a MSV of 0.7 or above. Only those GO terms exhibiting an FDR < 0.05 after Bonferroni Correction for multiple testing and a fold enrichment above 2 are shown. Table S10 (Associated to Figure S7A). List of genes constituting each GO term in Figure S5. GO term enrichment analysis of genes from clusters of mock-inoculated plants at the IN area with a MSV of 0.7 or above. Only those GO terms exhibiting an FDR < 0.05 after Bonferroni Correction for multiple testing and a fold enrichment above 2 are shown. Table S11 (Associated to Figure S7B). List of genes constituting each GO term in Figure S6. GO term enrichment analysis of genes from clusters of mock-inoculated plants at the OUT area with a MSV of 0.7 or above. Only those GO terms exhibiting an FDR < 0.05 after Bonferroni Correction for multiple testing and a fold enrichment above 2 are shown. Table S12. Primers used in this study and primer concentration for RT-qPCRs. Table S13. RT-qPCR results in numeric format along with Cp values of Targets and Cp value of Reference housekeeping gene. Dataset 1 Click here to access/download Dataset Table S1_Primers.xlsx Dataset 2 Click here to access/download Dataset Table S2_Statistics.xlsx Dataset 3 Click here to access/download Dataset Table S3_DEGs.xlsx Dataset 4 Click here to access/download Dataset Table S4_DEGsINvsOUT.xlsx Dataset 5 Click here to access/download Dataset Table S5_GOsINvsOUT.xlsx Dataset 6 Click here to access/download Dataset Table S6_GeneListClusters_IN.xlsx Dataset 7 Click here to access/download Dataset Table S7_GOsClusters_IN_PtoAvrRpm1.xlsx Dataset 8 Click here to access/download Dataset Table S8_GOsClusters_OUT_PtoAvrRpm1.xlsx Dataset 9 Click here to access/download Dataset Table S9_GOsDEGs.xlsx Dataset 10 Click here to access/download Dataset Table S10_GOsCluster_IN_Mock.xlsx Dataset 11 Click here to access/download Dataset Table S11_GOsCluster_OUT_Mock.xlsx Dataset 12 Click here to access/download Dataset Table S12_Primers.xlsx Dataset 13 Click here to access/download Dataset Table S13_qPCR_RawData.xlsx