Precision transcriptomics of viral foci reveals the spatial regulation of immune-signaling genes and identifies RBOHD as an important player in the incompatible interaction between potato virus Y and potato.

Whereas the activation of resistance (R) proteins has been intensively studied, the downstream signaling mechanisms leading to the restriction of the pathogen remain mostly unknown. We studied the immunity network response conditioned by the potato Ny-1 gene against potato virus Y. We analyzed the processes in the cell death zone and surrounding tissue on the biochemical and gene expression levels in order to reveal the spatiotemporal regulation of the immune response. We show that the transcriptional response in the cell death zone and surrounding tissue is dependent on salicylic acid (SA). For some genes the spatiotemporal regulation is completely lost in the SA-deficient line, whereas other genes show a different response, indicating multiple connections between hormonal signaling modules. The induction of NADPH oxidase RBOHD expression occurs specifically on the lesion border during the resistance response. In plants with silenced RBOHD, the functionality of the resistance response is perturbed and the spread of the virus is not arrested at the site of infection. RBOHD is required for the spatial accumulation of SA, and conversely RBOHD is under the transcriptional regulation of SA. Using spatially resolved RNA-seq, we also identified spatial regulation of an UDP-glucosyltransferase, another component in feedback activation of SA biosynthesis, thus deciphering a novel aspect of resistance signaling.

INTRODUCTION

Plants have evolved sophisticated mechanisms to perceive pathogen attack and effectively respond. Besides passive defenses, plants possess a multilayered actively induced immune system. The first most general layer, pathogen‐associated molecular‐pattern‐triggered immunity (PTI), is based on the recognition of pathogen components, common for a number of pathogens, by the plant pattern recognition receptors. The more specific layers of immunity are mediated by intracellular resistance (R) proteins (Jones and Dangl, 2006). R proteins confer the recognition of pathogen‐derived effectors and initiate effector‐triggered immunity (ETI). One of the manifestations of a successful immune response is hypersensitive response (HR)‐conferred resistance, where restriction of pathogens to the infection site is associated with a form of localized programmed cell death (PCD) (Künstler et al., 2016). HR‐conferred resistance is preceded by a series of biochemical and cellular signals. One of the earliest hallmarks of HR is the rapid and intense production of reactive oxygen species (ROS) (Balint‐Kurti Peter, 2019). Salicylic acid (SA) is required for the restriction of pathogens during HR in various pathosystems, including viruses (Mur et al., 2008; Künstler et al., 2016; Calil and Fontes, 2017) such as tobacco mosaic virus (TMV) (Chivasa et al., 1997; Chivasa and Carr, 1998) and potato virus Y (PVY) (Baebler et al., 2014; Lukan et al., 2018). The effectiveness of downstream events in ETI is also regulated by jasmonic acid (JA) and ethylene (ET) (Baebler et al., 2014). The PTI layer of immunity arising from plant–virus interaction is much less understood. It has been shown only recently that double‐stranded RNA (dsRNA) is a viral pathogen‐associated molecular pattern (PAMP) that triggers a signaling cascade similar to the PTI in other pathosystems (Niehl et al., 2016). Other hormones have also been shown to play important roles in plant immunity (Pieterse et al., 2012; Li et al., 2019). The activation of the immune signaling network results in the induced expression of actuators of defense, such as pathogenesis‐related protein 1 (PR‐1) and beta‐glucosidases, yet their function is not fully understood (Morant et al., 2008; Breen et al., 2017). HR‐associated PCD was shown to restrict pathogen spread in some (hemi)biotrophic pathosystems (Dickman and de Figueiredo, 2013; Lee et al., 2015b). It has been shown that it is not required for resistance in several viral pathosystems (reviewed in Künstler et al., 2016), however, including the potato–PVY interaction (Lukan et al., 2018).

Precise temporal and spatial coordination of induced signaling pathways is required to successfully restrict the pathogen with minimal damage to the host (reviewed in Künstler et al., 2016). The concentric spread, typical of many viruses, results in foci containing cells at different stages of infection (Yang et al., 2007; Rupar et al., 2015). The immune response signal is, however, transferred to the surrounding tissue, resulting in a gradient of response in the surrounding cells (Dorey et al., 1997; Havelda and Maule, 2000; Maule et al., 2002; Yang et al., 2007).

Potato virus Y (PVY), a member of the genus Potyvirus from the family Potyviridae, is the most harmful virus of cultivated potatoes (Karasev and Gray, 2013) and is among the most economically important viruses infecting plants (Quenouille et al., 2013). The most studied type of resistance to PVY is HR‐conferred resistance (Karasev and Gray, 2013). In potato cv. Rywal, HR is initiated by the Ny‐1 gene and manifests as the formation of necrotic lesions on inoculated leaves (Szajko et al., 2008). We have shown previously that the transcriptional dynamics of genes known to participate in the immune response is crucial for the efficient resistance response and that SA is the key component in the orchestration of these events (Baebler et al., 2014).

In order to dissect the mechanisms acting downstream of the activation of the Ny‐1 R protein in the resistance response, we analyzed the processes at work on the biochemical, ultrastructural and gene expression levels in the cell death zone and surrounding tissue in a spatiotemporal manner. To evaluate the position of SA in the signaling cascade, we analyzed in parallel the responses of SA‐depleted NahG plants.

RESULTS

The transcriptional response of immune signaling‐related genes is diversely spatially regulated

Previously, we performed a non‐targeted transcriptional analysis of inoculated homogenized leaves at 1, 3 and 6 days post‐inoculation (dpi) in potato cv. Rywal leaves infected with PVY (Baebler et al., 2014). The results of our several independent experiments showed that responses can vary between plants, which is also related to the number of lesions formed on the leaves. We thus hypothesized that a more detailed spatial analysis of responses would provide better insights into the processes involved in the resistance response downstream of Ny‐1 activation. In response to PVYN‐Wilga, the phenotypes of potato (cv. Rywal, hereafter NT) and a transgenic counterpart, depleted in the accumulation of SA (NahG‐Rywal, hereafter NahG, with SA reduced to <10% of native values; Baebler et al., 2014), were as previously described (Baebler et al., 2014; Lukan et al., 2018; Figure S1). Although resistance was successful in NT and the virus was restricted to the inoculated leaves, the depletion of SA rendered NahG plants susceptible (catechol independent; Baebler et al., 2014), allowing the virus to spread throughout the plant accompanied by the development of lesions (Figure S1).

We developed a protocol for the precise sampling of tissue sections surrounding the lesion at different stages of lesion development (early visible lesions at 3 dpi and fully developed lesions at 5–6 dpi; Figure 1a), suitable for transcriptomics and hormonomics analyses. As the quantities of RNA obtained from individual small tissue sections (approximatey 400 cells; Figure S2) are limited, we first selected 23 genes from our microarray data set according to their biological function and responsiveness (Baebler et al., 2014; Table S1), including genes participating in redox potential homeostasis (CAT1, PRX28, RBOHD, RBOHA, RBOHC, TRXO and TRXH), cell death (MC3 and LSD1), RNA silencing (AGO2 and SAHH), MAP kinase cascade (WIPK and MKP1), SA (TGA2), JA (13‐LOX, 9‐LOX, ACX3), ET (ERF1) or auxin (ARF2) signaling, actuators of defense (BGLU and PR1B) or a marker for starch metabolism (GBSSI) (Figure 1c). In parallel with the gene expression of those genes, the relative quantity of PVY RNA was measured by quantitative PCR (qPCR; Data S1 and S2).

Figure 1

Experimental set‐up.

(a) Sections containing lesion (section A) and surrounding tissues (sections B, C and D, consecutive 1‐mm strips distal to section A) were sampled at two different time points. Tissue surrounding early visible lesions was sampled in two perpendicular directions, here marked with 1 and 2, respectively. For electron microscopy, tissue sections with a lesion and surrounding cells were sampled. In mock‐inoculated plants, four adjacent tissue sections of the same size as used in inoculated leaves were excised from different plants. Scale bar: 1 mm.

(b) PVYN‐GFP accumulation around fully developed lesions in non‐transgenic (NT, top) or in SA‐depleted (NahG, bottom) plants. From left to right: PVYN‐GFP accumulation (green), chlorophyll fluorescence (red), overlay of chlorophyll fluorescence and PVYN‐GFP accumulation. Scale bar: 2 mm.

(c) Plant immune signaling‐related genes selected for the initial transcriptomic spatiotemporal response analysis. Potato virus Y (PVY) accumulation and expression of genes involved in reactive oxygen species (ROS), MAP kinase signaling, programmed cell death (PCD), salicylic acid (SA), jasmonic acid (JA), ethylene (ET), auxin and starch metabolism were analyzed (bold border). Signaling proteins are presented as orange, enzymes as red and transcription factors as blue ovals, respectively; metabolites are presented as diamonds and genes as squares. ACX3, acyl‐CoA oxidase; AGO2, argonaute 2; ARF2, auxin response factor 2; BGLU2, 1,3‐β‐glucosidase; CAT1, catalase 1; ERF1, potato ethylene responsive transcription factor 1a; GBSS1, granule‐bound starch synthase 1; HSP70, heat shock protein 70; 9‐LOX, 9‐lipoxygenase; 13‐LOX, 13‐lipoxygenase; LSD1, zinc‐finger protein LSD1; MC3, metacaspase 3; MKP1, mitogen‐activated protein kinase phosphatase 1; PR1B, pathogenesis‐related protein 1b; PRX28, peroxidase 28; RBOHA, potato respiratory burst oxidase homolog A; RBOHC, potato respiratory burst oxidase homolog C; RBOHD, potato respiratory burst oxidase homolog D; SAHH, adenosylhomocysteinase S‐adenosyl‐l‐homocysteine hydrolase; TGA2, bZIP transcription factor family protein; TRXH, thioredoxin H; TRXO, thioredoxin O; WIPK, potato wound‐inducible protein kinase. For gene IDs and primer information, see Table S3.

First, we analyzed fully developed lesions that developed 5–6 days following inoculation with PVYN‐Wilga (Figure 1a). The highest quantity of viral RNA was found in the center of the lesion in both potato genotypes, but the expression pattern of the genes changed with distance from the center of the lesion, and in some cases the response differed between the two genotypes (Figure 2). The response did not differ if the plant was inoculated with different viral strains, however (PVYNTN and PVYN‐GFP; Figure S3). As the lesion diameter is somewhat larger in NahG plants at the stage when we collected the fully developed lesions (up to twofold larger in NahG compared with NT plants; Figure S2; Data S3), note that some of the response measured in section B of the fully developed lesions of NahG plants is captured within section A in NT plants. We collected the samples as soon as the dark lesion was visible in order to minimize this effect. Such early sampling was also necessary to only follow the effect of the depletion of SA and not also the multiplication of PVY (Figures S2 and S3). The expression of HSP70 (see Figure 1c for full gene names), encoding a protein involved in the stability of R proteins (Kadota and Shirasu, 2012), is elevated in the center of viral foci (section A) compared with the adjacent sections (sections B, C and D) only in non‐transgenic (NT) plants (Figure 2). Both BGLU2 and PR1B are also under the regulation of SA (Figure 2). Although these genes do indeed show strong upregulation close to the center of viral infection (more than 10‐fold compared with the distal tissue, section D) during the resistance response, this response differs markedly in the SA‐depleted system. Whereas the spatial response of BGLU2 in NahG plants is completely lost, the peak of the PR1B response occurs further away from the foci of viral infection (reaching a peak in section C; Figure 2). The regulation of ERF1, involved in ET signaling, resembles that of PR1B in both genotypes. The spatial profile of starch synthase GBSS1, an indicator of primary metabolism state, shows that starch metabolism is downregulated in the center of the lesion, and its proximity, whereas in section D its expression status is closer to that of the mock‐inoculated leaf tissue (Figure 2, denoted with empty circles at the end of the x‐axis). This is not true in NahG plants, however, where the expression of GBSS1 is downregulated even further away from the lesion (Figure 2).

Figure 2

The lack of SA leads to the diverse spatial transcriptional regulation of immune‐related genes.

Transcriptional profiles of the selected genes were monitored in non‐transgenic (NT) and SA‐depleted (NahG) plants following inoculation with PVYN‐Wilga at the stage of fully developed lesions (5–6 dpi). Distance from the lesion center, marked as positions A, B, C and D, is plotted on the x‐axis and relative gene expression is plotted on the y‐axis. Expression values in mock‐inoculated tissue sections are shown as empty circles at the end of x‐axis. Asterisks denote the statistical significance (P < 0.05) of gene expression at positions B, C and D, compared with that at position A. Spatial profile models are shown as thick black lines, with 95% confidence interval bands presented in gray. Relative gene expression values within individual lesions are presented with colored symbols connected by a line. Full gene names are given in the legend of Figure 1.

We also show that the expression of genes involved in JA biosynthesis, in the maintenance of redox potential and in cell death is spatiotemporally regulated, which is presented in more detail in the next sections. For some genes the transcriptional response in the cell death zone and surrounding tissue depends upon SA, with a lack of SA leading either to no response (for BGLU2, ACX3, RBOHA, RBOHD, TRXH and GBSS) or to a changed spatiotemporal profile of the response (for ERF1, PR1B, CAT1 and PRX28) (Figures S4 and S5; Data S1). The rest of the genes analyzed are not SA‐dependent, showing either no spatial regulation at all or responding in an asynchronous manner between biological replicates (ARF2, RBOHC, 13‐LOX, WIPK, MKP1, AGO2, SAHH, LSD1, MC3, 9‐LOX, HSP70 and TGA2); Figures S4 and S5; Data S1), meaning that most probably their regulation is even more finely tuned and that more detailed analysis (e.g. with single‐cell resolution) would be required.

Next, we wanted to implement our approach in the earlier stages of viral infection. As the first robust time point for spatial analysis of responses, we identified the early visible lesions that develop at 3 dpi. These lesions are visible only when leaves are transilluminated (Figure 1a). Indeed, when comparing the response profiles of different lesions and profiles within the same lesion in perpendicular directions, the results were consistent and enabled a reliable biological interpretation (Figures S6 and S7; Data S2). Comparison of spatial response profiles in both time points enabled the identification of unexpected regulation. For example, the JA signaling module is not uniformly spatiotemporally regulated. The 13‐LOX gene, encoding 13 lipoxygenases at the start of the oxylipin branch, leading to JA biosynthesis, is induced in the center of viral foci at the stage of fully developed lesions (Figure S8). In agreement with the hypothesis that JA and SA are antagonistic, the response is enhanced in NahG plants. The expression pattern of the ACX3 gene, involved in later steps of JA biosynthesis (Figure 1c), is regulated similarly, except for the diminished regulation at the stage of fully developed lesions of NahG plants (Figure S8). This shows that besides JA signaling, SA signaling is also required for the spatiotemporal response of the ACX3 gene. Our results thus indicate that as yet unknown connections between different components of immunity signaling exist.

The ultrastructural features in the cell death zone and surrounding tissue are not dependent on SA

Our next hypothesis was that SA depletion influences processes that lead to changes in ultrastructural features in and surrounding the site of viral entry, thus rendering the immune response ineffective. In the center of lesions the cells were in the final stages of the cell death process, and only the cell wall and a small volume of electron‐dense material were visible (Figure S9). The transition between these and normal cells was limited to a narrow layer of cells of different types (Figure 3), both at the stage of early and fully developed lesions. In parenchymal and bundle sheath cells cytoplasmic vacuolization, an autophagy feature, was in progress (Figure 3a,c,f). Many chloroplasts were deformed or degraded (Figure 3b,f). No differences were observed between both the genotypes studied (Figure 3), however, indicating that the features of the processes observable at this level (PCD, autophagy, ER stress, etc.) do not depend on SA.

Figure 3

Ultrastructural features of the hypersensitive response are not dependent on salicylic acid. The same features were observed at the border of the cell death zone in both non‐transgenic (a–c) and NahG plants (d–f).

(a) Ameboidal nucleus and vacuolized cytoplasm in bundle sheath cell.

(b) Electron‐dense cytoplasm, disorganized thylakoid (distortion of grana with larger interthylakoidal spaces) in chloroplast and detachment of plasmalemma from the cell wall (white arrows) of the companion cell.

(c) Highly vacuolized cytoplasm and accumulation of electron‐dense material inside (white arrows) and outside (black arrow) the cell.

(d) Detachment of plasmalemma from the cell wall (white arrows) in bundle sheath cell.

(e) Degradation of chloroplasts and increased endoplasmic reticulum in mesophyll cell.

(f) Detachment of plasmalemma (white arrows), vacuolization of cytoplasm (black arrows), beginning of chromatin condensation and degradation of chloroplasts in bundle sheath cell. CH, chloroplast; CW, cell wall; CYT, cytoplasm; DC, destroyed cell; ER, endoplasmic reticulum; ES, extracellular space; GA, Golgi apparatus; M, mitochondrion; N, nucleus; P, peroxisome; PM, plasmalemma; V, vacuole.

The expression of genes involved in redox‐state maintenance is tightly regulated across the spatiotemporal scale studied

We next focused our research on genes involved in the generation and quenching of ROS as redox potential changes are one of the first responses to pathogen recognition, both in PTI and ETI (Figure 1c). We studied three RBOH genes involved in the generation of apoplastic ROS and leading to redox changes in the cytoplasm. On the quenching side, we analyzed gene expression of CAT1, located in the peroxisome, and apoplastic peroxidase PRX28. Additionally, we monitored the gene expression of two cytoplasmic redox potential sensor proteins that were implicated in the regulation of immune signaling, thioredoxins H and O (TRXH and TRXO). The three RBOH genes investigated each demonstrate distinct spatiotemporal regulation. RBOHC is not regulated in any of the genotypes studied (Figure 4). RBOHA is strongly induced in the center of viral foci already in the early stage of lesion development, whereas this activation is attenuated in NahG plants. RBOHD, however, responds in section A in the early stages, whereas later on it forms a peak of expression close to the border of the lesion itself in the majority of lesions (Figure 4). RBOHD is also under the regulation of SA as this gene is expressed close to the limit of quantification in NahG plants. PRX28 is induced earlier in NahG plants, and at the stage of fully developed lesion its expression peaks in section A in NT plants and in section B in NahG plants. CAT1 is induced only in NahG plants and is induced even further away from the center of viral foci compared with PRX28. At the stage of early lesions it peaks in section C, and at the stage of fully developed lesions it peaks even further from the center of the lesion, in section D. The expression of TRXO was below the limit of quantification, and therefore it was not possible to follow its response. Interestingly, however, TRXH responds strongly and early only in NT plants.

Figure 4

The induction of expression of RBOHD is limited to the cells immediately surrounding the infected tissue.

Spatial expression profiles of selected genes for non‐transgenic (NT) and SA‐depleted (NahG) plants after inoculation with potato virus Y at the stage of early visible lesions (3 dpi, left, averaged profile of the two perpendicular directions) and fully developed lesions (5–6 dpi, right). For the full gene names and full description of the graphs, see the legends to Figures 1 and 2 , respectively.

We further wanted to inspect how changes in transcriptional activity of ROS‐related genes are reflected in the accumulation of H2O2. We followed H2O2 production in relation to lesion size at several time points following virus inoculation. In contrast to the visible lesion growth pattern (Figure 5a), the H2O2 accumulation area around the lesions was significantly larger in NahG plants after 5 dpi (Figure 5b). Interestingly, the stained area continuously increased in NahG plants, whereas it reached a plateau in NT plants, indicating that ROS production is limited around the lesions in plants exhibiting resistance, but not in SA‐depleted plants (Figure 5c).

Figure 5

H2O2 accumulation is continuously greater during infection with potato virus Y in potato plants perturbed in SA signaling.

(a) Lesion diameter was significantly larger in the NahG plants from 5 dpi on compared with non‐transgenic (NT) plants.

(b) DAB staining indicating H2O2 activity (visible as the brown precipitate) was significantly wider in NahG plants from 5 dpi onwards, compared with NT plants.

(c) Difference in diameter of DAB staining and the diameter of the corresponding lesion, showing that in NahG plants the H2O2‐positive area outside the lesion is wider in comparison with NT plants. The representative areas of leaf images and comparison of lesions and DAB staining at 8 dpi are shown under the plots (scale bar: 10 mm). Asterisks represent statistically significant differences between the two genotypes on the same day (P < 0.05). Error bars represent standard errors.

Functional RBOHD contributes to the blocking of systemic PVY spread

RBOHD showed a specific spatiotemporal expression profile, peaking at the border of the virus amplification zone (Figure 4). To validate the involvement of ROS signaling and specifically the RBOHD gene in virus restriction, we constructed transgenic plants of cv. Rywal with silenced RBOHD gene (shRBOHD; Figure S10a). We inoculated two independent transgenic lines of shRBOHD (lines 13 and 14) with PVYN‐GFP as the least virulent PVY strain among those used in this study (Lacomme et al., 2017). RBOHD is, in selected transgenic lines, silenced to approximately 50% of its native transcriptional activity (Figure 6b). The virus spread systemically and lesions appeared in upper non‐inoculated leaves (Figures 6a and S10c,d). The time of first lesion appearance in upper non‐inoculated leaves was comparable with the time of first lesion appearance in NahG plants, but the viral RNA quantities were lower than those found in NahG plants (Figure S10c,d). Fully functional RBOHD thus contributes to the long‐distance movement of the virus.

Figure 6

The virus efficiently multiplies in RBOHD‐silenced plants.

(a) Lesions on an upper non‐inoculated leaf of RBOHD‐silenced plant of transgenic line 14 (right) and non‐transgenic counterpart (left) at 17 dpi. Scale bar: 1 cm.

(b) Relative RBOHD expression in inoculated leaves at 6 dpi.

(c) Relative potato virus Y (PVY) RNA abundance (relative to the lowest detected viral level, shown on a logarithmic scale) at 6 dpi in NT, NahG and shRBOHD plants. P, plant number; L, leaf number.

(d) Number of lesions per cm2 as they appeared on the youngest inoculated leaves of plants (n = 6–10) of different genotypes at a particular dpi. Error bars represent the standard errors of the number of lesions per cm2. Individual measurements are available in Data S4.

We also checked whether the virus can multiply to a larger extent in the inoculated leaves of shRBOHD plants by detecting viral RNA abundance and by counting the number of infection sites. The viral abundance negatively correlated (r = −0.85) with RBOHD expression, taking into account all of the genotypes analyzed (Figures 6b,c, S10b,c and S11a). Next, we followed the number of lesions appearing on the inoculated leaves of the RBOHD‐silenced transgenic lines after inoculation with PVY to study whether the reduction of RBOHD expression results in a higher number of PVY infection sites. We confirmed this hypothesis as the number of PVY foci on inoculated leaves of shRBOHD plants was higher compared with NT plants (Figures 6d and S11b–d; Data S4). Therefore, we concluded that the virus can infect cells more easily in shRBOHD plants and that it also multiplies to a larger extent in the inoculated leaves of these plants. Fully functional RBOHD thus contributes to blocking the initiation of viral infection, the local multiplication of the virus and the long‐distance movement of the virus.

RBOHD is involved in the tight spatial regulation of SA accumulation during the resistance response

The ROS were shown to be involved both upstream and downstream of SA signaling (Herrera‐Vasquez et al., 2015). To determine whether RBOHD‐generated ROS are involved in the regulation of SA signaling in the resistance response, we analyzed SA induction in shRBOHD plants. In the pool of lesions (section A), SA was induced to the same extent as in RBOHD‐silenced plants (approximately 10‐fold compared with mock samples; Figure 7a; Data S5), whereas in the surrounding cells (section B) the accumulation was lower in non‐transgenic plants but not in RBOHD‐silenced plants. As the concentration of SA in lesions on transgenic plants did not differ significantly from the concentration in lesions on non‐transgenic plants (Figure 7a), we conclude that RBOHD‐generated ROS are not involved in the induction of SA biosynthesis in potato‐PVY interaction. The tight spatial regulation of SA biosynthesis was lost in shRBOHD plants, however. Additionally, we know that RBOHD gene activity is under the regulation of SA as the expression of this gene is much lower in NahG plants after virus inoculation (Figures 6b and S10b; Data S1).

Figure 7

Functional RBOHD contributes to the spatial regulation of salicylic acid accumulation in potato resistance against potato virus Y.

(a) Salicylic acid (SA) content was measured in the early stage of lesion development in non‐transgenic (NT) and RBOHD‐silenced (shRBOHD) plants within the lesion (section A) and in the tissue surrounding the lesion (section B) in comparison with mock‐inoculated sections. Pools of 50–100 tissue sections from between six and eight plants were analyzed. Detailed data are shown in Data S5.

(b) Expression of UDP‐glucosyltransferase (UGT76D1) within the lesion (section A), its close surroundings (section B) and sections of mock‐inoculated plants of NT and shRBOHD plants, as determined by RNA‐seq analysis. The bar graph presents the average normalized read counts (CPM) of the UGT76D1 gene across three biological replicates. Asterisks represent statistically significant differences among comparisons (false discovery rate‐adjusted P < 0.05).

(c) SA synthesis and RBOHD are interconnected through a regulatory feedback loop. SA regulates RBOHD expression and, on the other hand, RBOHD is, through UGT76D1, involved in the spatial distribution of SA accumulation. SA regulates the chloroplastic production of reactive oxygen species. Signaling proteins are presented as orange ovals and metabolites as red diamonds. JA/ET, jasmonic acid/ethylene signaling; Ny‐1, R protein; ROS, reactive oxygen species.

We have further explored the mechanisms underlying spatial regulation of SA accumulation on the border of the lesion. We performed RNA‐seq analysis of the tissue within the lesion (section A) and the adjacent tissue (section B) at the time of early visible lesions for all three genotypes: NT, NahG and shRBOHD. This analysis provided insights into processes regulated by SA or RBOHD (Figures 8 and S12; Data S6). The majority of processes that are differentially regulated between sections A and B overlap between NahG plants and shRBOHD plants (Figure 8). This was expected, as the expression of RBOHD is repressed in NahG plants (Figure 6b). We also demonstrate that some of the genes and even processes/protein groups are specifically regulated only when the RBOHD gene is silenced, however. Most notable is the regulation of MYB transcription factors (three regulated only in the absence of RBOHD, and not in NT or NahG transgenic plants; Data S6) and the regulation of protein degradation (the specific regulation of several proteases as well as ubiquitin SCF and RING E3 ligases; Data S6). This shows that RBOHD has some additional roles that are independent of SA.

Figure 8

RBOHD contributes to the spatial regulation of processes in potato resistance signaling against viruses.

Differentially regulated processes between the lesion (section A) and the surrounding tissue (section B) of non‐transgenic plants (NT), plants with perturbed SA accumulation (NahG) or plants with reduced RBOHD expression (shRBOHD). Only statistically significant (false discovery rate‐corrected Q < 0.05) enriched gene sets determined by gene set enrichment analysis (GSEA) in at least one comparison are presented. Process groups and individual processes together with mapman bin names and the numbers of genes in the respective bin are shown. ‘+’, induced processes; ‘−’, repressed processes; blanks denote that a current process was not statistically significantly enriched.

We further inspected the RNA‐seq results for any evidence of a feedback loop indicating the role of RBOHD in the spatial regulation of SA accumulation. We have identified a UDP‐glucosyltransferase (UGT76D1, Sotub10g024000) that is strongly spatially deregulated in RBOHD‐silenced plants and could be directly involved in the process (Figures 7b and S13). In line with the evidence that RBOHD is repressed in NahG plants, a similar expression pattern is also observed in these plants. In Arabidopsis, UGT76D1 is implied in the feedback activation loop of SA biosynthesis (Huang et al., 2018). Overexpression of this gene led to the accumulation of SA. Here, we similarly show that the silencing of RBOHD results in the induction of the UGT76D1 gene in section B, which may underlie the higher levels of SA detected in section B of transgenic lines (Figure 7b).

DISCUSSION

The mechanism of the activation of R proteins and their interactions with effectors has been the subject of intensive research in recent years (Baggs et al., 2017; Kourelis and Van Der Hoorn, 2018; Zhou and Zhang, 2020). To work towards sustainable resistance in the field, where the plant is also exposed to other environmental stressors, understanding the downstream processes leading to pathogen arrest is of importance to balance the trade‐offs between growth and immunity, however. To date, most of the studies of those processes have been performed in the model plant Arabidopsis (Piquerez et al., 2014). Although some findings can be transferred to crop species using orthology (Lee et al., 2015a; Ramšak et al., 2018), this is not always the case, as in the recently reported redundancy of PAD4 in Solanaceae, for example (Gantner et al., 2019). Therefore, it is important to also perform the studies in crop plants, such as the potato.

Pathogen‐infected leaf tissue comprises a heterogeneous mixture of host cells in different stages of defense response as a result of pathogen and defense signal movement from the primary infected cell (Maule et al., 2002). Recently, spatially resolved transcriptome profiling in plant tissues has gained importance in unraveling complex regulatory networks (Giacomello et al., 2017; Shulse et al., 2019). The progressive and asynchronous effects of viral infection on host gene expression in relation to the spatial distribution of the virus have been investigated previously, but in compatible interactions (Aranda et al., 1996; Escaler et al., 2000; Maule et al., 2002; Yang et al., 2007). In resistance, however, the spatial distribution of a limited number of components has been studied only in model plants (Antoniw and White, 1986; Mur et al., 1997; Betsuyaku et al., 2018; Giolai et al., 2019). Our approach to the analysis of small tissue sections presents a step forward in studying the resistance response against the virus, and enabled the identification of a novel key player and the interconnectivity of components of immune signaling.

In the immune response, ROS are generated in two phases (Herrera‐Vasquez et al., 2015). The first, low‐amplitude and transitory phase occurs within minutes of infection. The ROS generated at this stage are mostly apoplastic, tightly linked to the post‐translational activation of plasma membrane RBOH NADPH oxidases (e.g. AtRBOHD and AtRBOHF) and cell wall peroxidases, and are attributed to PTI. (Torres, 2010; Shapiguzov et al., 2012; Kadota et al., 2019). The second, high‐amplitude and sustained phase takes place a few hours after infection and depends on ROS generation in multiple compartments, including the apoplast, chloroplasts, mitochondria and peroxisomes. This second phase is linked to ETI and requires the transcriptional activation of RBOH genes (Shapiguzov et al., 2012). The role of RBOH proteins as the principal generator of ROS after pathogen attack was mostly studied in compatible interactions with bacterial and oomycete pathogens, using mutant and transgenic lines (Allan et al., 2001; Torres et al., 2002; Yoshioka et al., 2003; Peer et al., 2011), but also in symbiotic interaction with rhizobia (Yu et al., 2018). Our spatiotemporal analysis of responses revealed that potato RBOHD (ortholog of AtRBOHD) is transcriptionally induced at the border region between virus‐replicating and healthy cells around the cell death zone in the resistance response (Figure 4). Several studies have pointed to the essential role of ROS also in HR‐conferred virus resistance (reviewed in Hernández et al. 2016), including in one of the most studied viral pathosystems between Nicotiana tabacum (tobacco) and the tobacco mosaic virus (TMV) (Mur et al., 1997; Liao et al., 2015). Although the induction of RBOHD in the potato–PVY interaction has been shown previously (Otulak‐Kozieł et al., 2019), its role in the resistance response has not been studied functionally. We directly confirmed its essential role in potato resistance signaling, as RBOHD‐silenced transgenic plants were not able to arrest the spread of the virus (Figure 6a,c). Besides the effect on cell‐to‐cell movement, both SA and RBOHD also have an effect on the efficiency of viral infection. The number of established foci of viral infection was the highest in NahG plants followed by RBOHD‐silenced plants, whereas the lowest number of viral foci appeared in non‐transgenic plants (Figure 6d). These results could imply a stronger impact of SA in immune signaling compared with RBOHD; however, one must note that NahG reduces the levels of SA to <10% of native values, whereas RBOHD is silenced to approximately 50% of its native transcriptional activity in our transgenic lines (Figure 6b).

Although the oxidative burst following pathogen recognition occurs in the apoplast, being generated mostly by RBOH proteins and cell wall peroxidases, pathogen‐induced ROS can be also produced in other cellular compartments like mitochondria, chloroplasts and peroxisomes (reviewed in Sharma et al., 2012; Mignolet‐Spruyt et al., 2016). Our comparative analysis of NT and NahG genotypes shows that the extensive production of ROS corresponds to the spread of the virus and/or PCD (Figure 5; Lukan et al., 2018). The extensive ROS response in later stages of infection is most likely of non‐apoplastic origin, as ROS production did not correspond to the expression of RBOH genes in NahG plants (Figures 5, 6and5, 6). This is also supported by the expression pattern of ROS quenchers CAT1 and PRX28, which are induced more distantly from the virus foci in NahG compared with NT plants (Figure 4). As SA was shown to be required for the generation of mitochondrial ROS during pathogen infection (Liao et al., 2015), chloroplasts are its most probable source in our system. Transcriptional downregulation of photosynthesis genes leads to the generation of chloroplast ROS, which were proposed as a signal orchestrating PCD (Zurbriggen et al., 2010; Su et al., 2018). In our pathosystem the transcriptional downregulation of photosynthesis genes was observed (Baebler et al., 2014), and thus we can assume that, consequently, chloroplastic ROS are generated. Our results of H2O2 staining (Figure 5a,b) are in agreement with the hypothesis that chloroplastic ROS stimulates localized cell death. Moreover, as the stained area was limited around the lesions in plants exhibiting resistance, but not in SA‐depleted plants (Figure 5c), we suggest that SA regulates chloroplastic H2O2 production.

The studies in Arabidopsis could not unambiguously determine the role of SA in the initiation of pathogen‐induced PCD. In some pathosystems PCD can be regulated by SA, but not in all (Huysmans et al., 2017; Radojičić et al., 2018). Here, we determined that features of HR PCD are not dependent on SA in our studied pathosystem (Figures 3, 6and 6d). Moreover, this study provides the first detailed spatiotemporal ultrastructural analysis of HR PCD in potato following Potyvirus infection. We observed typical features of PCD and autophagy on the border of early and fully developed lesions (Figure 3). Similar features were previously observed in different pathosystems resulting either in susceptible or resistant plant–Potyvirus interaction (Hinrichs‐Berger et al., 1999; Otulak and Grazyna, 2012; Zielińska et al., 2012; Choi et al., 2016), supporting the hypothesis that the mechanisms of PCD are independent of the resistance mechanism (Künstler et al., 2016). This is also in line with a recent study revealing that the activation of NLR proteins can directly lead to cell lysis through the insertion of resistosome into plasma membrane (Wang et al., 2019).

Reactive oxygen species (ROS) were also proposed to be a central component of a self‐amplifying loop that regulates the interaction balance between different phytohormones such as SA, JA and ET (reviewed in Torres, 2010). SA regulates RBOHD‐dependent ROS production in Arabidopsis (reviewed in Liu and He, 2016). The same is true for our pathosystem, as the induction of RBOHD expression was significantly reduced after PVY inoculation in NahG plants (Figures 6b and S10b; Data S1). In Arabidopsis, RBOHD knockout mutant plants accumulated higher levels of SA following interaction with pathogens (Pogany et al., 2009). On the other hand, Chaouch et al. (2012) did not detect any difference in SA accumulation in the RBOHD mutant compared with wild‐type Arabidopsis (Chaouch et al., 2012). We showed that in the potato–PVY interaction, SA biosynthesis is not controlled by RBOHD‐generated ROS, as the concentration of SA in transgenic plants with silenced RBOHD did not differ significantly from SA concentration in NT plants at the site of viral foci. RBOHD is, however, involved in the spatial control of its accumulation (Figure 7a,c), confirming the existence of complex regulatory feedback loops. Indeed, UDP‐glucosyltransferase (UGT76D1) was identified as a component of the feedback activation loop of SA biosynthesis (Huang et al., 2018) and our results show that it is spatially repressed by RBOHD signaling, explaining the regulation of SA accumulation (Figure 7b,c).

For efficient immunity the plant needs to block pathogen multiplication and spread, and at the same time limit the extent of damaged tissue and energy consumption. Thus the tight spatial regulation of these processes is of utmost importance. It has been long assumed that positive regulators act at the HR site and negative regulators act in the surrounding areas, but the molecular evidence for this premise is mostly lacking through a lack of functional zonation studies (Huysmans et al., 2017). We show here that RBOHD and SA regulate the immune response in a spatiodiverse manner, leading to the fine‐tuning of both spatial and temporal responses.

EXPERIMENTAL PROCEDURES

Plant material

Potato (Solanum tuberosum L.) plants cv. Rywal and NahG‐Rywal (Baebler et al., 2014) and derived transgenic lines (see below) were grown and inoculated with PVYNTN (isolate NIB‐NTN, AJ585342), PVYN‐Wilga (PVYN–Wi; EF558545) or PVYN605, tagged with green fluorescent protein (PVYN‐GFP) (Jakab et al., 1997; Dietrich and Maiss, 2003; Rupar et al., 2015) or mock inoculated, as described by Baebler et al. (2009).

qPCR gene expression analysis in tissue sections

For gene expression analyses lesions at two different stages were sampled (see Figure 1a). For each lesion, four tissue sections were sampled (Figure 1a). Tissue sections were stored in 100 µl of RNAlater RNA Stabilization Solution (ThermoFisher Scientific, https://www.thermofisher.com). Altogether, five independent experiments were performed, each comprising of approximately 20 plants per experimental group. For each experimental group, between five and 10 lesions from different plants were sampled (Table S2).

Standard gene expression analysis procedures were optimized for the analysis of small leaf sections. Fixed tissue sections were homogenized using Tissue Lyser (Qiagen, https://www.qiagen.com), followed by RNA isolation using the RNeasy Plant Micro Kit (Qiagen), DNAse treatment, quality control and reverse transcription.

The expression of 23 genes involved in different steps of immune signaling (Figure 1c) was analyzed and normalized to the expression of two validated reference genes, COX and 18S, as described previously (Petek et al., 2014; for primer and probe information and full experimental details, see Table S3), to eliminate technical variability, including the effect of different quantities of RNA in tissue sections. The standard curve method was used for relative gene expression quantification using quantgenius (http://quantgenius.nib.si; Baebler et al., 2017).

Statistical modeling of tissue section qPCR gene expression data

Prior to a statistical analysis of the gene expression dataset for tissue sections, data from lesions without viral amplification detected in section A (i.e. necrosis cause by mechanical inoculation) were filtered out. Relative gene expression values were next standardized to the 97th quantile of the maximum value across all samples (for inoculated and mock samples in both studied genotypes). The multiple linear regression method, as implemented in the r package stats (R Core Team, 2020), was used to fit a quadratic polynomial model, with genotype and position (tissue section) or measured distance from position A (initial viral foci tissue surrounding) included as factors. An analysis of variance (ANOVA) was additionally used to analyze the effects of position and genotype on gene expression. For each gene analyzed, spatial expression profiles were determined and plotted with 95% confidence intervals. For data for the early visible lesions, calculations were performed for both perpendicular directions and their average, respectively. Expression changes between genotypes at different times post‐inoculation were further visualized using coefficient plots, with linear coefficients plotted on the x‐axis and quadratic coefficients plotted on the y‐axis. The significance of comparisons between the levels of the factors was obtained through the contrast method as implemented in the r package limma.

Gene expression in whole leaves

For the gene expression analysis in whole leaves of different genotypes, including newly constructed transgenic lines (see the experimental set‐ups below), whole leaves were sampled. RNA was isolated and reverse transcription and qPCR was performed as described previously (Baebler et al., 2014), with normalization to COX only. Where applicable, a Student’s t‐test (in excel) was used to compare treatments.

RNA‐seq analysis in tissue sections

For RNA‐seq, early visible lesions were sampled from PVYN‐Wilga‐inoculated leaves of cv. Rywal, NahG‐Rywal and shRBOHD (line 14) plants. Thirty tissue sections of the lesion (A) and its immediate surroundings (B; Figure 1a) were collected and pooled separately (pool A and pool B). Sections of mock‐inoculated plants were prepared as described above and pooled (with approximately 300 sections per pool). The pools were stored and homogenized as described above. RNA extraction and library preparation are described in Method S1. RNA‐seq was performed on the HiSeq platform (Illumina, https://www.illumina.com) using 150‐bp paired‐end reads at Novogene (https://en.novogene.com) or LC Sciences (https://www.lcsciences.com). Quality control was performed in clc genomics workbench 12.0 (Qiagen). The merging of overlapping pairs, the mapping of reads to the potato genome (Petek et al., 2020) and read counting were performed in both clc genomics and star. Differential expression analysis was performed in r (R Core Team, 2020; version 3.2.2), using the r package limma (Ritchie et al., 2015), trimmed mean of M values (TMM) normalization and voom function. Adjusted P values below 0.05 were considered statistically significant. For full experimental details, see Method S1. Raw and normalized RNA‐seq data were deposited to the Gene Expression Omnibus (GEO, accession number GSE142002, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE142002).

Gene set enrichment analysis (GSEA; Subramanian et al., 2005) was performed to search for groups of genes involved in the same processes that were significantly (false discovery rate‐corrected Q < 0.05) altered by virus inoculation, using MapMan ontology as the source of the gene sets (obtained from GoMapMan database; Ramšak et al., 2014).

Construction of short hairpin RNA transgenics

Primers were designed based on a consensus sequence obtained by sequence alignment of sequences of RBOHD genes (Sotub06g025550 and Sotub06g025580) (Figure S14). cDNA obtained from cv. Rywal was used as a template for the amplification of the 436‐bp‐long fragment of the RBOHD gene. The product was cloned into the pENTR D‐TOPO plasmid and sequenced. The fragment was further transferred to pH7GWIWG2, sequenced and electroporated into Agrobacterium tumefaciens LBA4404 that was used for the transformation of cv. Rywal plantlets (for full experimental details, see the legend to Figure S14). Silencing efficacy was analyzed in one leaf per transgenic line by qPCR using shRBOHD primers (Table S3) as a target, as described above.

Phenotypization

Lesions that developed on virus‐inoculated leaves of cv. Rywal, NahG‐Rywal and two short hairpin transgenic lines with silenced RBOHD (lines 13 and 14) inoculated with PVYN‐GFP or mock inoculated were counted from 3 to 14 dpi (for the numbers of plants tested for each experiment, see Data S4). The number of lesions was calculated per cm2 of leaf area. For the same plants, the times of lesion appearance on non‐inoculated leaves was also recorded.

Viral quantity was analyzed in the same genotypes as above at 6 dpi in the inoculated leaves and at 24 and 39 dpi in non‐inoculated leaves (using between one and three leaves from up to five plants). Sample preparation and qPCR analysis were performed as described above using PVY as a target gene. The experiment was repeated three times.

The lesion diameter was measured in the images of individual leaves of cv. Rywal and NahG‐Rywal, inoculated with PVYN‐Wilga or PVYNTN, in photoshop cs3 (Adobe, https://www.adobe.com ) and standardized to leaf diameter. The experiment was repeated twice.

H2O2 staining

Leaves (three leaves from three plants per genotype/treatment) of Rywal and NahG‐Rywal plants, inoculated with PVYNTN or mock inoculated, were sampled 3–8 days after inoculation and stained with 3,3′‐diaminobenzidine (DAB) staining solution (for full experimental details, see Method S2). Diameters of lesions and DAB staining areas were measured on spatially calibrated images in fiji (Schindelin et al., 2012). The difference in the diameter of lesion and the DAB area was measured for each individual lesion. Data were analyzed in r 3.5.0 (R Core Team, 2020), and statistical significance was tested with anova and Tukey’s post‐hoc test (P = 0.05).

Transmission electron microscopy

Ultrastructural changes following PVYNTN and PVYN‐Wilga inoculation of the cells in and around the early visible and fully developed lesions were observed in Rywal and NahG‐Rywal using transmission electron microscopy (TEM; Philips CM 100). A tissue section including a lesion and surrounding cells was cut from inoculated leaves and prepared as described in Method S3. To estimate the lesion size in each genotype, cross‐sections of the leaves were analyzed by TEM. The average lesion size was determined by counting the number of dead cells in the diameter of approximately 15 lesions for each genotype.

Hormonal measurements

Both SA and JA contents were determined in non‐transgenic cv. Rywal and two short hairpin transgenic lines with silenced RBOHD (lines 13 and 14) at 4 dpi after PVYN‐Wilga or mock inoculation. Between 60 and 130 sections containing initial virus foci (section A) and surrounding tissue (section B) were sampled from the third inoculated leaf from between six and eight plants for each genotype. Hormones were isolated and measured by gas chromatography coupled with mass spectrometry (GC‐MS) as described by Križnik et al. (2017).

Funding information

This research was financially supported by the Slovenian Research Agency (research core funding no. P4‐0165 and projects J4‐7636, J4‐1777 and 1000‐15‐0105) and grant no. 2017/27/B/N29/01803 from the Polish National Science Centre.

CONFLICT OF INTEREST

The authors declare that they have no conflicts of interest.

AUTHOR CONTRIBUTIONS

ŠB, MPN, KG, TL AND JH designed the research. ŠB, MPN, KG, TL, MTŽ, AK, KS, BD, AC, SP, KM and MK performed the research. TL, ŠB, AK, MTŽ, AB, MZ, KG, MK and SP analyzed the data. TL, ŠB, AK, MZ, KG, MK, MPN, MTŽ, BD, AC and JH contributed to the writing or revision of the article.

Open research badges

This article has earned an Open Data badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at {https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE142002}.

Figure S1. Lesions appear systemically if Ny‐gene‐conferred resistance is perturbed in SA signaling.

Figure S2. Lesion diameter is twofold larger in NahG compared with NT plants following inoculation with PVY.

Figure S3. Transcriptional response of selected immune signaling‐related genes and relative viral RNA abundance are similar after inoculation with different viral strains.

Figure S4. Spatiotemporal gene expression of immune signaling genes.

Figure S5. 3D representation of spatial gene expression response to PVY inoculation at the time of fully developed lesions.

Figure S6. 3D representation of spatial gene expression response to PVY inoculation at the time of early visible lesions.

Figure S7. 3D animated spatial gene expression of immune signaling genes around early visible lesions.

Figure S8. 9‐LOX and 13‐LOX branches of the oxylipin pathway respond differentialy in potato resistance against PVY.

Figure S9. Programmed cell death in the center of a fully developed lesion after PVY infection.

Figure S10. The virus is present in symptomatic non‐inoculated leaves of plants of RBOHD‐silenced transgenic lines.

Figure S11. Performance of RBOHD‐silenced plant inoculated leaves after PVY inoculation?

Figure S12. Spatial transcriptional response is similar in shRBOHD and NahG genotypes.

Figure S13. Spatial expression of UDP‐glucosyltransferase (UGT76D1) in different genotypes.

Figure S14. Construction of short hairpin RNA transgenics.

Table S1. Gene expression changes of selected genes in whole‐leaf homogenates at different time points after inoculation with PVY.

Table S2. Independent experiments performed for targeted spatiotemporal gene expression analyses.

Table S3. Selected genes and corresponding primers and probe sequences for expression analyses with quantitative PCR.

Methods S1. RNA‐seq analysis.

Methods S2. H2O2 staining and imaging.

Methods S3. Transmission electron microscopy.

Data S1. Standardized qPCR gene expression data.

Data S2. Results of statistical modeling of tissue section gene expression data.

Data S3. Individual lesion sizes.

Data S4. Number of lesions on inoculated leaves of different potato genotypes following PVY inoculation.

Data S5. Hormonal contents in tissue sections of different potato genotypes.

Data S6. Results of spatial RNA‐seq analysis of Rywal, NahG and shRBOHD genotypes following PVYN‐Wilga inoculation.

Click here for additional data file.