Transcriptome Analysis Shows Activation of Stress and Defense Responses by Silencing of Chlorophyll Biosynthetic Enzyme CHLI in Transgenic Tobacco

In the present study, we have shown the transcriptional changes in a chlorosis model transgenic tobacco plant, i-amiCHLI, in which an artificial micro RNA is expressed in a chemically inducible manner to silence the expression of CHLI genes encoding a subunit of a chlorophyll biosynthetic enzyme. Comparison to the inducer-treated and untreated control non-transformants and untreated i-amiCHLI revealed that 3568 and 3582 genes were up- and down-regulated, respectively, in the inducer-treated i-amiCHLI plants. Gene Ontology enrichment analysis of these differentially expressed genes indicated the upregulation of the genes related to innate immune responses, and cell death pathways, and the downregulation of genes for photosynthesis, plastid organization, and primary and secondary metabolic pathways in the inducer-treated i-amiCHLI plants. The cell death in the chlorotic tissues with a preceding H2O2 production was observed in the inducer-treated i-amiCHLI plants, confirming the activation of the immune response. The involvement of activated innate immune response in the chlorosis development was supported by the comparative expression analysis between the two transgenic chlorosis model systems, i-amiCHLI and i-hpHSP90C, in which nuclear genes encoding different chloroplast proteins were similarly silenced.


Introduction
Viruses induce varying symptoms in their hosts by affecting the morphological and physiological alterations of the hosts. Chlorosis is one of the most common symptoms in plant-virus interaction. Reduced productivity along with significant yield loss are common features of chlorotic plants. Virus-induced chlorosis, represented by the distortion and dysfunction of chloroplasts, results from the reduced chlorophyll production and photosynthetic gene expression [1][2][3][4]. Although the major strategies for protecting plants from virus diseases comprise the use of resistant cultivars and the control of vector organisms, regulation of symptom development could contribute to alleviating crop loss. Therefore, understanding of precise mechanisms underlying chloroplast dysfunction and

RNA Sequencing, Mapping, and Identification of Differentially Expressed Genes (DEGs)
Three-week-old i-amiCHLI transgenic lines and non-transformed tobacco (SR1) plants were treated with a dexamethasone (Dex) solution or a control solution. When they were cultured for an additional week, the Dex-treated CHLI transgenic lines showed chlorosis ( Figure 1A) as described previously [19]. RNA extracted at 24 h post-treatment was subjected to RNA-seq, which gave around 23 M reads/sample. From 74 to 90% of those reads of 12 samples were mapped to the tobacco reference transcriptome (Supplementary Table S1). The overall similarity of transcriptomes in all samples was examined using a principal component analysis (PCA). All three Dex-treated C-1 samples are readily distinguishable from the control-treated C-1 samples, and Dex-and control-treated non-transformed SR1 samples ( Figure 1B). However, control-treated C-1 plants have also shown an apparent difference from the non-transformed SR1 plants ( Figure 1B). This difference between non-transformant SR1 and control-treated i-amiCHLI plants could be due to a slight downregulation of the CHLI gene in the control-treated plants ( Figure 1C). To verify the reduced CHLI expression in the control-treated C-1 plants, we compared the expression levels by qRT-PCR of CHLI and representative genes that show clear up-and down-regulation found in the RNA-seq analysis (see below), encoding the nuclear-protein systemic acquired resistance deficient 1 (SARD1) and chloroplastic light-harvesting chlorophyll a/b-binding protein (LHC a/b), respectively, using primers shown in Supplementary Table S2. Unexpectedly, the control-treated C-1 plants showed comparable expression levels of CHLI, SARD1, and LHC a/b to SR1 plants, while Dex-treatment induced reproducible CHLI silencing accompanied by the up-and down-regulation of SARD1 and LHC a/b, respectively ( Figure 1D,E). The SR1 plants did not respond to Dex-treatment in the expression of the genes examined. These results indicate that the different transcriptome profiles of control-treated C-1 plants from those in SR1 plants cannot be attributed to the difference in the temporal CHLI expression levels.
The gene expression of Dex-treated i-amiCHLI plants (group CD; CD1, CD3, and CD4) was compared in a pairwise manner using DESeq2 with all three control groups, control-treated line C-1 (group CC; CC1, CC2, and CC4), Dex-treated SR1 (group SD; SD1, SD3, and SD4), and control-treated SR1 (group SC; SC2, SC3, and SC4), in which we have never observed any chlorosis in numbers of experiments. The criteria for DEGs selection were as follows: the absolute values of log2 (FC) larger than 1 with the adjusted p-values lesser than 0.05. The MA plots suggest similar numbers of genes were differentially expressed in the three comparisons (Supplementary Figure S1). The analyses of differentially expressed mRNA transcripts annotated with the Arabidopsis AGI codes in three different comparisons above, CD vs. CC, CD vs. SD, and CD vs. SC, identified 1481, 2088, and 1425 up-regulated genes, respectively, and 1046, 2237, and 937 down-regulated genes, respectively (Figure 2A,B, Supplementary Tables S3 and S4).
Because the control-treated C-1 plants have shown some difference in the overall expression pattern from both control SR1 plants (group SC and SD) and Dex-treated C-1 plants (group CD), we thought that the use of commonly up-or down-regulated genes in three comparisons for gene ontology (GO) analysis was inadequate. Indeed, among DEGs detected in CD vs. SD and CD vs. SC comparisons, 438 upregulated and 381 downregulated genes were not detected in CD vs. CC comparison (Figure 2A,B). It is possible that genes essential for chlorosis development would be up-or down-regulated in both CD and CC groups. Therefore, we used in GO enrichment analysis 1428 (45.9%) and 982 (31.8%) genes, which were commonly up-and down-regulated in at least two comparisons, respectively (Figure 2A,B, Supplementary Tables S3 and S4). In the CC vs. SC comparison, we found 1592 upregulated and 2206 downregulated genes (Supplementary Tables S3 and S4), which were also subjected to GO enrichment analysis (see below). Table S4). It is notable that 276 genes in the profile 27 were up-regulated in CC group but downregulated in CD as compared to SR1 plants (SC and SD) ( Figure 2D and Supplementary  Table S4). The STEM analysis enabled us to pick up more DEGs than those we found in the pairwise analyses, in which some important genes in chlorosis development could have been masked by the unexpectedly affected transcriptome of control-treated C-1 plants.  , i-amiCHLI-1 (C-1) and i-amiCHLI-2 (C-2) plants were grown for three weeks, sprayed with 0.01% Tween-20 containing 0.5% ethanol (control) or 50 µM Dex (Dex), and photographed at 7 days post-treatment (dpt). Scale bars indicate 1 cm. (B) Principle component analysis (PCA) of RNA-seq data. RNA preparations from 1-dpt control-treated SR1 (SC2, SC3, and SC4) and Dex-treated SR1 (SD1, SD3, and SD4), control-treated C-1 (CC1, CC2, and CC4), and Dex-treated C-1 (CD1, CD3, and CD4) were analyzed using DESeq2. (C) The relative expression level of CHLI in all RNA-seq samples from DESeq2 analysis. (D,E) Expression analysis of SARD1 (D) and LHC a/b (E). The expression levels of these genes relative to a fixed standard sample (square) were quantified in eight individual plants each of control-(circles) and Dex-treated (triangles) i-amiCHLI-1 plants, and in four individual plants each of control-(crosses) and Dex-treated (diamonds) SR1 plants. The relative expression values were plotted with those of CHLI. and Dex-treated C-1 (CD1, CD3, and CD4) were analyzed using DESeq2. (C) The relative expression level of CHLI in all RNA-seq samples from DESeq2 analysis. (D,E) Expression analysis of SARD1 (D) and LHC a/b (E). The expression levels of these genes relative to a fixed standard sample (square) were quantified in eight individual plants each of control-(circles) and Dex-treated (triangles) i-amiCHLI-1 plants, and in four individual plants each of control-(crosses) and Dex-treated (diamonds) SR1 plants. The relative expression values were plotted with those of CHLI.  In addition to DESeq2, we employed the Short Time-series Expression Miner (STEM) program [21] for the analysis of four groups of RNA-seq data. Although STEM was designed to analyze RNA-seq data from time-course samples, we thought it would be useful for our analysis because it can cluster genes based on the expression profiles. We subjected the relative normalized count data from four plant/treatment groups, SC, SD, CC, and CD, of 15,498 genes annotated with Arabidopsis AGI codes to the analysis. The STEM clustering algorithm grouped all the genes into 50 different profiles (Supplementary Figure S2), and among them, significantly enriched profiles (p-values lesser than 0.05) were highlighted with different color for each cluster of profiles. The profiles 18,40,42,29, and, 30 were selected for the up-regulated genes in Dex-treated C-1 plants ( Figure 2C), containing a total 3568 up-regulated DEGs (Supplementary Table S3). Profiles 27,9,11,26,34,23, and, 31 were selected for the downregulated genes ( Figure 2D,E), containing a total of 3582 downregulated DEGs (Supplementary Table S4). It is notable that 276 genes in the profile 27 were up-regulated in CC group but downregulated in CD as compared to SR1 plants (SC and SD) ( Figure 2D and Supplementary Table S4). The STEM analysis enabled us to pick up more DEGs than those we found in the pairwise analyses, in which some important genes in chlorosis development could have been masked by the unexpectedly affected transcriptome of control-treated C-1 plants.

GO Enrichment Analysis of the DEGs
We performed GO enrichment analysis with the above chlorosis-related DEGs: 1428 up-and 982 down-regulated DEGs from the pairwise analysis (Supplementary Table S5) and 3635 up-and 3583 down-regulated DEGs from STEM analysis (Supplementary Table S6). GO terms enriched by DEGs detected by two different methods were compared, indicating that about half of the GO terms are common to the both. Table 1 shows the selected GO terms for biological processes enriched by the lists of up-regulated DEGs. The GO terms enriched by the up-regulated DEGs comprise innate immune response, defense response, response to wounding, response to hormones [salicylic acid (SA), jasmonic acid (JA), and abscisic acid (ABA)], and plant-type hypersensitive response (Supplementary Figure S3) were present (Table 1 and Supplementary Tables S5 and S6). Moreover, SA biosynthetic genes were also found to be up-regulated (Table 1). These results suggest that the reduction of CHLI levels in chloroplast induces defense response, including the activation of some phytohormone pathways, which can be regarded as candidates of chlorosis-accelerating genes.
The GO terms (biological process) enriched by the downregulated genes include photosynthesis, photoinhibition, pigment metabolic process, and plastid organization ( Table 2, Supplementary Tables S5  and S6). Several GO terms that correspond to some primary metabolism genes were also downregulated. These GO terms include glycine metabolic process, carbohydrate metabolic process, lipid metabolic process, and vitamin metabolic process. Some cellular process genes annotated with the GO terms of cellular homeostasis, chaperone-mediated protein folding, and cellular component organization were downregulated (Table 2, Supplementary Tables S5 and S6), which is consistent with the upregulation of cell death-related genes ( Table 1, Supplementary Tables S5 and S6). Also, genes involved in cellular response to chemical stress, response to oxidative stress, and response to heat were downregulated ( Table 2). As mentioned before, we compared the transcriptomes in CC plants that did not show any visible chlorosis with those of SC plants and found significant transcriptomic differences between them. Among them, 954 out of 1592 upregulated genes and 1336 out of 2206 downregulated genes were also detected in more than one of CD vs. CC, CD vs. SD, and CD vs. SC comparisons (Supplementary Figure S4A,B). Comparison of GO enrichment analyses with the chlorosis-related DEGs and those in CC plants revealed that about half of GO terms were common to these analyses (Supplementary Figure S4C,D). Interestingly, GO terms related to SA and cell death pathways were enriched by chlorosis-related DEGs but not by DEGs of CC plants, while both sets of DEGs enriched GO terms associated with JA, ABA, and plant immunity (Supplementary Table S3).

Cell Death Assay in i-amiCHLI Plants
Interestingly, the GO term of cell death was significantly enriched in the up-regulated DEGs. As shown in Figure 3A, in Dex-treated i-amiCHLI plants the cell death-related genes have shown clear upregulation compared to the control-treated i-amiCHLI and non-transformant SR1 plants. This finding encouraged us to detect cell death in i-amiCHLI lines. Cell death in upper and lower leaves was determined separately since the lower leaves had shown more severe chlorosis. Micrographs of the trypan blue-stained leaf tissue prove cell death in lower leaves of the Dex-treated C-1 and C-2 plants, albeit with the lesser extent and less shrunken morphology ( Figure 3B: h,l), compared to those in the positive control, in which tomato bushy stunt virus (TBSV) P19 had transiently been expressed for two days ( Figure 3B: m). These observations suggest that the cell death events in Dex-treated i-amiCHLI plants are sporadic and slower than that of hypersensitive response, although the genes annotated with GO term of plant-type hypersensitive response were upregulated therein ( Table 1). The upper leaves of Dex-treated C-1 plants showed cell death of a lesser extent and those of Dex-treated C-2 plants did not show any staining ( Figure 3B: k), which was similar to the cases in control-treated i-amiCHLI and non-transformant SR1 plants ( Figure 3B: a-f,i,j). These observations suggest that the cell death observed in this study is an age-dependent event. Electrolyte leakage assay supported the notion above: an age-dependent, sporadic, and slow cell death ( Figure 3C).

Detection of ROS Production in i-amiCHLI Plants
The dysfunction of chloroplasts often results in the production of ROS. Especially, when the regulation of chlorophyll biogenesis is impaired, the intermediate compounds serve as strong photosensitizers that lead to the accumulation of ROS [7,12]. Our GO enrichment analysis failed to detect upregulation of genes related to ROS response but rather some of them were downregulated (Table 2). However, we found subsets of ROS response-related genes were significantly up-regulated in one of the Dex-treated C-1 plants ( Figure 4A). Also, the other two Dex-treated C-1 plants were differentially clustered from control-treated C-1 plants and SR1 plants in the clustering analysis of ROS response-related gene expression of all the plants ( Figure 4A). These results prompted us to detect the H 2 O 2 production as a representative of ROS. In the positive control samples where intense light was applied, the accumulation of brownish precipitate in chloroplasts was found by diaminobenzidine (DAB) staining ( Figure 4B: m). Intense DAB staining of chloroplasts was barely detected in Dex-and control-treated SR1 plants, and control-treated C-1 and C-2 plants ( Figure 4B: a-f,i,j). Clear DAB staining of chloroplasts was detected in leaves of Dex-treated C-1 and C-2 plants at 1-dpt ( Figure 4B: g,h,k,l), unlike the samples from 7-dpt when cell death was observed (data not shown). These results suggest that the silencing of CHLI induces ROS production in the chloroplasts.

Comparison of Chlorosis-Related Transcriptome Changes in i-amiCHLI and i-hpHSP90C Transgenic Plants
Previously, we described the transcriptome changes preceding the chlorosis developments and relevant biological events in another chlorosis model, i-hpHSP90C transgenic plant system [20], which is similar to the present i-amiCHLI plants but differs in the silencing target. To find out the common features in the chlorosis development in these two transgenic systems, we compared the results from GO enrichment analyses with up-and down-regulated genes. We found 105 and 179 common GO terms enriched in the analyses with up-and down-regulated genes, respectively ( Figure 5A,B), suggesting the involvement of common pathways in the chlorosis development in two chlorosis model systems. The common GO terms of upregulated genes include defense, cell death, SA, JA, ABA, and stress responses ( Figure 5C and Supplementary Figure S3), and those of down-regulated genes include chloroplast biogenesis, photosynthesis, and development ( Figure 5D). The results suggest that similar pathways including the defense responses accompanying cell death have roles in chlorosis development in two chlorosis models with different chlorosis triggers.

Detection of ROS Production in i-amiCHLI Plants
The dysfunction of chloroplasts often results in the production of ROS. Especially, when the regulation of chlorophyll biogenesis is impaired, the intermediate compounds serve as strong photosensitizers that lead to the accumulation of ROS [7,12]. Our GO enrichment analysis failed to light was applied, the accumulation of brownish precipitate in chloroplasts was found by diaminobenzidine (DAB) staining ( Figure 4B: m). Intense DAB staining of chloroplasts was barely detected in Dex-and control-treated SR1 plants, and control-treated C-1 and C-2 plants ( Figure 4B: a-f,i,j). Clear DAB staining of chloroplasts was detected in leaves of Dex-treated C-1 and C-2 plants at 1-dpt ( Figure 4B: g,h,k,l), unlike the samples from 7-dpt when cell death was observed (data not shown). These results suggest that the silencing of CHLI induces ROS production in the chloroplasts.  In addition to those common GO terms, we also found some GO terms enriched specifically in either of the two systems. To confirm the specific nature of those up-and down-regulated genes, we further analyzed them by clustering (Supplementary Figure S5). The results showed that the CD4 sample exhibited similar expression patterns to HD2 and HD3, which are analyzed as the chlorosis-developing samples in our previous study (Supplementary Figure S5C-L). Assuming that chlorosis develops faster in the HSP90C system, the results of the comparative analysis suggest that CD4 plants had preceded CD1 and CD3 on the way to develop chlorosis, and that CD1 and CD3 had not undergone some of the processes leading to the chlorosis development. The results suggest that the expression of genes assigned to those GO terms (Supplementary Figure S5C-L) does not essentially differ between the two model systems, and thus, support the idea that chlorosis development in two systems shares molecular mechanisms. By contrast, we also found that genes assigned to GO terms of "ERAD (endoplasmic reticulum-associated degradation) pathways" (Supplementary Figure S5A) and "Response to endoplasmic reticulum (ER) stress" (Supplementary Figure S4B) were specifically up-regulated in the HSP90C system: CD group samples including CD4 were differentially clustered from HD2 and HD3. To test the hypothesis that the silencing of CHLI leads to the downregulation of HSP90C, which drives transcriptome changes leading to chlorosis development, we examined the expression levels of six HSP90 family proteins in the RNA-seq data. In contrast to HD2 and HD3 samples, which showed downregulation of the silencing target HSP90C (HSP90.5) and upregulation of ER-localizing HSP90.7, CD group plants showed marginal decrease of HSP90C expression but no significant change in HSP90.7 expression (Supplementary Figure S6). Instead, CC and CD group plants showed reduced expression of cytoplasmic HSP90.1 (Supplementary Figure S6). However, it remains unknown whether or not such differences in HSP90 family proteins reflect the difference of molecular mechanisms underlying the chlorosis development. Further studies are necessary to elucidate the involvement of the ER-localizing and cytoplasmic HSP90 proteins in chlorosis development.

Comparison of Chlorosis-Related Transcriptome Changes in i-amiCHLI and i-hpHSP90C Transgenic Plants
Previously, we described the transcriptome changes preceding the chlorosis developments and relevant biological events in another chlorosis model, i-hpHSP90C transgenic plant system [20], which is similar to the present i-amiCHLI plants but differs in the silencing target. To find out the common features in the chlorosis development in these two transgenic systems, we compared the results from GO enrichment analyses with up-and down-regulated genes. We found 105 and 179 common GO terms enriched in the analyses with up-and down-regulated genes, respectively ( Figure  5A,B), suggesting the involvement of common pathways in the chlorosis development in two chlorosis model systems. The common GO terms of upregulated genes include defense, cell death, SA, JA, ABA, and stress responses ( Figure 5C and Supplementary Figure S3), and those of downregulated genes include chloroplast biogenesis, photosynthesis, and development ( Figure 5D). The results suggest that similar pathways including the defense responses accompanying cell death have roles in chlorosis development in two chlorosis models with different chlorosis triggers. In addition to those common GO terms, we also found some GO terms enriched specifically in either of the two systems. To confirm the specific nature of those up-and down-regulated genes, we further analyzed them by clustering (Supplementary Figure S5). The results showed that the CD4

Discussion
We have previously generated two inducible silencing systems for CHLI and HSP90C in tobacco and found that the silencing of these genes could stimulate the plants to develop chlorosis [18][19][20]. As discussed previously [20], these systems have advantages over the experimental systems with the virus-or viroid-infected plants because, in these transgenic lines, it is possible to analyze those plant cells that have committed to developing but not exhibited chlorosis. Utilizing the advantages of these systems, in our previous and current studies, an RNA-seq method was implemented to explore the early molecular changes, which lead to chlorosis development. Unlike previous studies analyzing the transcriptome changes in virus-and viroid-infected plants [3,22,23], we successfully identified transcriptome changes that precede the detectable chlorosis.
In the RNA-seq analysis, we found that the control-treated C-1 plants (CC group) were clustered differently from the control-and Dex-treated non-transformant SR1 plants (SC and SD groups) ( Figure 1B). Because the analysis of normalized count data in RNA-seq and qRT-PCR did not demonstrate the reduced CHLI expression in control-treated C-1 plants compared to SR1 plants, we have not identified the cause of altered transcriptome in control-treated C-1 plants. The qRT-PCR analysis of CHLI, SARD1, and LHC a/b indicated that the altered transcriptome of control-treated C-1 plants was most unlikely to result from reduced CHLI expression. The control-treated C-1 plants might have suffered some stress by the transgene during their growth. Nonetheless, the qRT-PCR analyses supported the qualitative difference between the transcriptomes of CC and CD groups, as up-and down-regulation of SARD1 and LHC a/b, respectively, were consistently observed only in Dex-treated C-1 plants ( Figure 1F,G).
Despite that the altered transcriptome of CC group plants increased the difficulty of RNA-seq data analysis, we identified several characteristic transcriptome changes in Dex-treated C-1 plants (CD group) through two types of analysis. One was pairwise analysis using DESeq2 as we used in the previous study [20], and the other is expression profile analysis using the STEM program. Because the latter could statistically be less stringent than the former, we carefully compared the GO terms enriched by the DEGs identified in those analyses. Like in the previous study of HSP90C-silenced plants [20], the downregulated were the genes involved in photosynthesis, different primary and secondary metabolic pathways, various cellular processes that control cellular homeostasis, and plant and cell growth, and upregulated were the genes involved in the triggering of innate immune response accompanying cell death, and also activation of different phytohormone signaling pathways. Also, we confirmed the ROS production and sporadic cell death in chlorosis-developing, Dex-treated i-amiCHLI transgenic plant lines, as we did in the Dex-treated i-hpHSP90C transgenic plant lines [20]. Therefore, the present study supports the idea from the previous study that (i) downregulation of chloroplast and photosynthesis-related genes (CPRGs) by retrograde signaling serves as the secondary cue for the chlorosis development, which follows the primary cue of the silencing of chloroplast protein genes with pivotal roles such as HSP90C or CHLI, and that (ii) upregulation of defense and cell death genes accelerate the chlorosis development through cell death induction. The comparative analysis of transcriptome changes in two models strongly supported this notion ( Figure 5).
The CPRGs downregulation has been extensively reported as a characteristic feature during the plant-virus interactions [3,[24][25][26]. Because CPRGs downregulation is quite usual in the chlorotic tissues, our current and previous studies [18][19][20]27,28] have consistently supported that the reduced expression of CPRGs preceding visible chlorosis is the primary pathway for chlorosis development. The present findings support the involvement of retrograde signaling (RS) in chlorosis induction induced by the CHLI silencing. The downregulation of the genes involved in maintaining cellular homeostasis, and plant and cell growth could be attributed to the RS activation, which reprograms transcriptome from growth and differentiation state to stress response state [29]. The RNA-seq data have supported the RS activation as the expression of transcription factors (TFs) involved in RS-mediated transcriptome changes was increased (Supplementary Table S8) [30]. We detected the H 2 O 2 production in chloroplasts on the next day of CHLI silencing induction. The upregulation of WRKY33 and WRKY40 suggests that ROS is a primary signal for the activation of RS pathways in the present chlorosis model. This notion is supported by a report that GUN4 and GUN5/CHLH generate singlet oxygen when bound with protoporphyrin IX [31], which is likely to accumulate in CHLI-silenced plants. Further studies would elucidate the mode of activation of the RS pathways in the present model system, and the involvement of RS in chlorosis development.
By upregulated genes, GO terms related to the innate immune response, cell death, and phytohormone signaling pathways were commonly enriched in two chlorosis model systems. In both systems, we confirmed the sporadic cell death. The mode of defense response activation remains to be studied, but phytohormone pathways could have roles. Particularly, genes involved in the regulation of SA biosynthetic pathways were found to be upregulated in both model systems. It is well-known that SA induces cell death in plants as a defense response [32,33]. In addition to cell death, we also observed ROS production that precedes visible chlorosis in both systems. Together, these results suggest that the induction of cell death is a common feature of plants with chloroplast dysfunction. Such cell death could be first triggered by the ROS production and then activated through the SA-ROS self-amplification loop [34]. Although the significance of SA and ROS in the cell death in our chlorosis model systems needs to be studied further, the comparison of CC vs. SC plants suggested the crucial role of SA and cell death in the development of visible chlorosis (Supplementary Table S3). As we discussed in the previous report, the cell death induction has barely been described in compatible interaction between plants and viruses [20]. Actually, CMV harboring Y-sat induces bright yellow symptoms [15,16] unlike the symptom-like phenotypes we observed in our systems. The difference in symptoms could result from the inhibition of cell death by a pathogens' counter-defense. In case of CMV, a viral factor involved in the cell death inhibition would be 2b protein, which reportedly have roles in the modulation of both SA-mediated resistance [35] and symptoms development [36].
In the comparative analysis of RNA-seq data from two chlorosis models, we also found that GO terms of "ERAD pathways" and "response to ER stress" were specifically enriched in i-hpHSP90C plants (Supplementary Figure S5A,B). A retrograde signal molecule, methylerythritol cyclodiphosphate, and SA reportedly activate the unfolded protein response (UPR) [37][38][39][40]. However, the role of UPR pathway in chlorosis development remains to be studied. Further study would provide us with an insight into the role of UPR in chlorosis development in the i-hpHSP90C model system, which may reveal the difference of two chlorosis model systems in the signaling pathways activating the major cues for chlorosis development: the downregulation of CPRGs and the activation of cell death. It is important to analyze the roles of ROS, SA, other phytohormones, and cell death event in the chlorosis development in the current model systems and others to understand the diverse mechanisms underlying the chlorosis. The outcome of such studies would lead us to identify genes that positively and negatively regulate the development of chlorosis, which could be novel targets for molecular breeding for crop protection.

RNA Isolation and Sequencing
We extracted total RNA using Sepasol-RNA 1 Super G (Nacalai Tesque, Kyoto, Japan) and treated the RNA with RNase-free DNase (Takara Bio, Kusatsu, Japan)as described previously [27]. RNA was extracted from six individual plants each of the Dex-and control-treated, SR1 and C-1 transgenic plants grown for 24 h after the Dex-or control-treatment. Then the RNA samples were analyzed using 2100 Bioanalyzer with RNA-nano chip (Agilent) to select three samples each with higher RNA quality from plant/treatment groups. TruSeq RNA Sample Preparation v2 kit (Illumina, Tokyo, Japan) was used for library construction. Sequencing of 100-bp single-end was performed at Nodai Genome Center using Hiseq2500 (Illumina).

Analysis of RNA-seq Data
RNA-seq data were acquired and trimmed as described previously [20] using bcl2fastq2 (Illumina) for acquisition, fastq_quality_trimmer for trimming. The transcript abundance estimation and expression value computation were done on the local Galaxy platform [41] using Salmon [42] with a reference transcriptome (Ntab-TN90-AYMY-SS_NGS.mrna.annot.fna downloaded from https://solgenomics.net/organism/Nicotiana_tabacum/genome) [43]. A previously described Transcript ID-AGI code table [20] was used to have the gene expression values with AGI codes, and 121,268 transcripts (64.02%) out of 189,413 transcripts of tobacco reference transcriptome were given AGI codes and included in the analysis [20]. For the identification of the DEGs, we used DESeq2, which provides a differential expression analysis method using negative binomial generalized linear models [44,45]. The "Gene Quantification" files from the Salmon analysis were subjected to DESeq2 analysis to give regularized log transformation of raw read count data. In the DESeq2 package, the "plotPCA" function is implemented, and this function was used to generate the PCA plots to visualize the sample-to-sample distances. The differentially expressed genes were selected under the following criteria: corrected p-value [P(adj)] less than 0.05 and the Log2(FC) values above 1.0 or below -1.0. The STEM software was used according to the method described by Liu and colleges [21] to find the detailed gene expression patterns of all the DEGs and relative normalized count data from DESeq2 were used for STEM analysis. GO enrichment analyses, preparation of heatmaps, and Venn diagrams were performed as described [20]. GOs with FDR lesser than 0.05 were regarded to be significant. The comparison of expression levels of selected genes was done using the relative normalized count data obtained in DESeq2 analysis.

Quantitative RT-PCR
Quantitative RT-PCR was conducted as described previously [20] using the primers in Supplementary  Table S2. Briefly, we extracted total RNA using the ISOSPIN Plant RNA Kit (Nippon Gene, Tokyo, Japan) from eight (transgenic) or four (SR1) biological replicates. Then, we synthesized cDNA using the M-MLV RTase (New England Biolabs, Tokyo, Japan). Real-time PCR was performed in the StepOnePlus Real-Time PCR system (Applied Biosystems, Thermo Fisher Scientific, Tokyo, Japan) using KAPA SYBR FAST qPCR master mix (Kapa Biosystems, via Nippon Genetics, Tokyo, Japan). The qPCR conditions were as follows: initial holding at 95 • C for 20 s, 40 cycles of 95 • C for 3 s, 60 • C for 30 s, followed by 95 • C for 15 s, 60 • C for 1 min, 95 • C for 15 s. The reaction specificity was examined by the melting curves. Each RNA preparation was analyzed with three technical replicates for each target and the reference gene. The relative expression levels of the genes of interest were calculated by the ∆∆CT method with EF1α as the internal reference.

Determination of Cell Death
For the determination of cell death, the trypan blue assay was conducted essentially as described [46] with a slight modification [20]. Briefly, three-week-old transgenic and control plants were Dex-or control-treated, grown for seven days, and observed for the phenotypic changes. Tobacco leaves with TBSV P19 [47] transiently expressed for 2 days served as a positive control for cell death [48]. Leaf disks were excised from three independent plants each of plant/treatment groups, stained and observed under a light microscope as described [20]. The cell death was quantitatively evaluated by measuring electrolyte leakage as described previously [49]. The cell death was evaluated with the relative conductivity in three triplicate experiments, as previously described [20]. We performed statistical analyses using SPSS (Version 17, IBM, Chicago, IL, USA) and Microsoft Office Excel 2016 (Microsoft Corporation, Albuquerque, NM, USA).

DAB (3,3'-Diaminobenzidine) Staining for Hydrogen Peroxide Detection
In situ detection of hydrogen peroxide (H 2 O 2 ), were performed by using the 3,3'-diaminobenzidine (DAB) staining as described previously [50] with slight modifications [20]. Briefly, three 6-mm leaf disks were excised from lower and upper leaves of three individual transgenic or control plants treated with Dex or control solution for 24 h, vacuum-infiltrated with 1 mg/mL DAB solution containing 0.05% Tween-20, incubated under light (70-100 µmolem −2 s −1 for 30 min for test plants and 250-300 µmolem −2 s −1 for 60 min for positive control plants) and additional 3.5 h in the dark on wet filter papers in Petri dishes, and de-colorized by boiling in 99.5% ethanol for 5-10 min and observed.