Transcriptional reprogramming caused by the geminivirus Tomato yellow leaf curl virus in local or systemic infections in Nicotiana benthamiana

Abstract

Background

Viruses have evolved to create a cellular environment permissive for viral replication in susceptible hosts. Possibly both enabling and resulting from these virus-triggered changes, infected hosts undergo a dramatic transcriptional reprogramming, the analysis of which can shed light on the molecular processes underlying the outcome of virus-host interactions. The study of the transcriptional changes triggered by the plant DNA viruses geminiviruses is potentially hampered by the low representation of infected cells in the total population, a situation that becomes extreme in those cases, like that of Tomato yellow leaf curl virus (TYLCV), in which the virus is restricted to phloem companion cells.

Results

In order to gain insight into how different the transcriptional landscapes of TYLCV-infected cells or whole tissues of TYLCV-infected plants might be, here we compare the transcriptional changes in leaf patches infected with TYLCV by agroinfiltration or in systemic leaves of TYLCV-infected plants in Nicotiana benthamiana. Our results show that, in agreement with previous works, infection by TYLCV induces a dramatic transcriptional reprogramming; the detected changes, however, are not equivalent in local and systemic infections, with a much larger number of genes differentially expressed locally, and some genes responding in an opposite manner. Interestingly, a transcriptional repression of the auxin signalling pathway and a transcriptional activation of the ethylene signalling pathway were detected in both local and systemically infected samples. A transcriptional activation of defence was also detectable in both cases. Comparison with the transcriptional changes induced by systemic infection by the geminivirus Tobacco curly shoot virus (TbSV) shows common subsets of up- and down-regulated genes similarly affected by both viral species, unveiling a common transcriptional repression of terpenoid biosynthesis, a process also suppressed by the geminivirus Tomato yellow leaf curl China virus.

Conclusions

Taken together, the results presented here not only offer insight into the transcriptional changes derived from the infection by TYLCV in N. benthamiana, but also demonstrate that the resolution provided by local and systemic infection approaches largely differs, highlighting the urge to come up with a better system to gain an accurate view of the molecular and physiological changes caused by the viral invasion.

Background

As intracellular parasites, viruses have evolved to create a cellular environment permissive for viral replication in susceptible hosts; for this purpose, viruses induce a re-wiring of the host’s physiology and development concomitant to the establishment of a successful infection. In plants, these virus-induced changes can be easily visualized and quantified, and infection by different viruses frequently produces some of a common array of symptoms, including stunting, chlorosis, and leaf curling. Possibly both enabling and resulting from these virus-triggered changes, infected hosts undergo a dramatic transcriptional reprogramming; the analysis of the modifications in the transcriptional landscape of the host upon the viral infection can shed light on the molecular processes underlying the outcome of virus-host interactions. Such transcriptional studies have proliferated in the past decade, possibly owing to technical advances allowing for in-depth sequencing, availability of genomic information, and the increased affordability of these approaches.

Geminiviruses are insect-transmitted DNA viruses causing severe diseases in crops worldwide, and currently pose a serious threat to food security; however, our understanding of the molecular basis of the infection is still partial, which limits the development of effective anti-geminiviral strategies for crop protection. The transcriptional changes triggered by the infection by geminiviruses have been studied in a number of plant-virus interactions [1,2,3,4,5,6,7,8]. When comparing the results obtained in these studies, few commonalities arise: plant hormone signaling pathways, especially those for jasmonates (JA) and brassinosteroids (BR), frequently appear as altered, although the direction of the change is not consistent [4, 5, 8, 9]; and cell cycle-related genes, which need to be reactivated in the virus-infected terminally differentiated cells to allow for viral DNA replication, are detected as differentially expressed in a couple of cases only [2, 9]. In an attempt to unveil the molecular basis of tolerance or recovery, Chen et al. (2013) and Gongora-Castillo et al. (2012) [3, 4] compared the transcriptome of susceptible or tolerant tomato cultivars infected with Tomato yellow leaf curl virus (TYLCV) and that of recovered and symptomatic leaves of pepper infected with Pepper golden mosaic virus (PepGMV), respectively, by RNA-seq; however, and somewhat surprisingly, only limited differences were detected in both cases.

One factor frequently neglected in these transcriptional studies of geminivirus-infected plants is the low representation of infected cells in the total population: in an average infection, only some cells will be supporting viral replication at a given time. This situation becomes extreme in those cases, like that of TYLCV, in which the virus is restricted to phloem companion cells. The global transcriptome of infected plants will be the average of that in all cells, infected and non-infected, and all tissues: this not only creates a serious dilution issue, but also averages the potential transcriptional responses in infected and systemic, uninfected cells, which could be opposite, hence generating potentially misleading results difficult to interpret. However, and since to date no transcriptional profile of isolated infected cells is available, the extent to which this could be different to that obtained from complete aerial organs of an infected plant is unclear.

In order to gain insight into how different the transcriptional landscapes of TYLCV-infected cells or whole tissues of TYLCV-infected plants might be, here we compare the transcriptional changes in leaf patches infected with TYLCV by agroinfiltration or in systemic leaves of TYLCV-infected plants in Nicotiana benthamiana by RNA-seq. Strikingly, our results show that, as expected and in agreement with previous works, infection by TYLCV induces a dramatic transcriptional reprogramming; the detected changes, however, are not equivalent in local and systemic infections, with a much larger number of genes differentially expressed locally, and some genes responding in an opposite manner in local and systemic samples. Interestingly, a transcriptional repression of the auxin signaling pathway and a transcriptional activation of ethylene signaling and defence responses were detected both in local and systemic infections. Despite more limited changes detected in the systemically infected samples, comparison with the transcriptional changes induced by systemic infection by the geminivirus Tobacco curly shoot virus (TbSV) [5] unveiled common subsets of up- and down-regulated genes similarly affected by both viral species. Among the common biological processes potentially affected by the transcriptional changes we find terpenoid biosynthesis as transcriptionally repressed; notably, another geminivirus species, Tomato yellow leaf curl China virus (TYLCCN), has been shown to suppress terpenoid biosynthesis [10], making it tempting to speculate that depletion of terpenoids might be a requirement for geminiviruses to establish a successful infection in nature. Taken together, the results presented here not only shed light on the transcriptional changes derived from the infection by TYLCV in N. benthamiana, but also demonstrate that the resolution provided by local and systemic infection approaches largely differs, highlighting the urge to come up with a better system to gain an accurate view of the molecular and physiological changes caused by the viral invasion.

Results

Transcriptional changes upon local infection by TYLCV in N. benthamiana

From the observation that only a fraction of cells support active replication by a geminivirus at a given time [11, 12] logically follows the idea that an accurate study of the cellular changes triggered by the viral invasion will require the isolation of the infected cells specifically, and their comparison with similar cells from an uninfected sample. Analysis of whole organs of systemically infected plants, which is the common practice due to the lack of a more precise approach, would presumably result in potential dilution and masking issues. In order to test this idea, we decided to compare the transcriptional changes detectable by RNA sequencing (RNA-seq) upon local or systemic infection in the model plant N. benthamiana by the geminivirus TYLCV.

The leaf patch infection results in the viral replication in most cells [11], which the virus must effectively manipulate to generate a permissive environment, hence serving as a good surrogate system to study the viral infection. For the local infection, we performed agroinfection in leaf patches of four-week-old N. benthamiana leaves and took samples at 6 days post-infiltration (dpi) (Fig. 1a), when the virus is still actively replicating. In order to exclude the potential effect of the bacteria on plant transcription, N. benthamiana leaves inoculated with an Agrobacterium clone containing the empty vector were used as control; three independent biological replicates were used in each case. RNA-seq was performed by Illumina sequencing as indicated in the methods section. The raw HiSeq reads were filtered and trimmed, and between 34 and > 51 million clean pair-end reads were obtained per sample; these clean reads were mapped to the N. benthamiana draft genome (v1.0.1) from the Sol Genomics Network (ftp://ftp.solgenomics.net/genomes/Nicotiana_benthamiana/assemblies/)) with a mapping rate between 89 and > 98% (Additional file 6 and Table 1). The PCA analysis of the three biological replicates for TYLCV and control samples is shown in Additional file 1: Figure S1.

Table 1 Summary of the RNA-seq results

Fig. 1

Transcriptional changes upon local infection by TYLCV in N. benthamiana. a Schematic representation of the experimental design. b Differentially expressed genes in the local infection are shown by volcano plot. The x-axis shows the log2 transformed gene expression fold change between infected and control samples. The y-axis indicates the negative log10 transformed adjusted p-values (FDR) of the differential expression test calculated by R package edgeR. The up-regulated and down-regulated genes are represented by red and blue dots, respectively. Pie chart shows the number of up−/down-regulated genes. c Validation of selected DEGs by qPCR. Values are the average of three biological replicates, relative to mock. NbACT was used as the normalizer. d Mapping of viral reads to the TYLCV genome. Three biological replicates (inner concentric circles) are represented; the upper side of each circle represents the virion (+) strand; the lower side of each circle represents the complementary (−) strand. ORFs are depicted in blue. Please note that the accumulated reads in the area containing the CP and V2 ORFs have been trimmed; the numbers of total reads are shown in Additional file 7: Table S2. e Validation of the expression of Rep and CP by qPCR. Expression values are relative to NbACT

A total of 7561 and 4289 down- and up-regulated genes, respectively, were identified in these locally infected samples (Fig. 1b; Additional file 6: Table S1). The RNA-seq results were validated by qPCR analysis of selected genes (Fig. 1c).

The clean reads were also mapped to the TYLCV genome, with an average of 88,550 reads per million (RPM) (Fig. 1d; Additional file 7: Table S2). Unexpectedly, reads for both strands of the virus were detected throughout the viral genome with uneven and non-perfectly symmetrical distribution, and not restricted to the described open reading frames (ORFs); the number of reads was much higher in the region of the genome containing the V2 and CP (late) genes (Additional file 7: Table S2). The accumulation of viral reads was confirmed by qPCR to detect expression of the viral genes encoding the Rep (Replication-associated protein) and the CP (capsid protein), contained in the complementary and the virion strand of the viral genome, respectively (Fig. 1e).

Transcriptional changes upon systemic infection by TYLCV in N. benthamiana

In order to analyze the transcriptional changes detectable during the systemic infection by TYLCV, we performed agroinfection of two-week-old N. benthamiana plants as described in [12], and took samples at 14 days post-infection (dpi) (Fig. 2a), when the virus is actively replicating in the apical leaves and the first symptoms have already appeared (Additional file 2: Figure S2). Apical leaves of N. benthamiana plants inoculated with an Agrobacterium clone containing the empty vector were used as control; three independent biological replicates were used in each case. RNA-seq was performed by Illumina sequencing as indicated in the methods section. The raw HiSeq reads were filtered and trimmed, and between 36 and > 52 million clean pair-end reads were obtained per sample; these clean reads were mapped to the N. benthamiana draft genome (v1.0.1) from the Sol Genomics Network (ftp://ftp.solgenomics.net/genomes/Nicotiana_benthamiana/assemblies/) with a mapping rate between 97 and > 98% (Additional file 6: Table S1). The PCA analysis of the three biological replicates for TYLCV and control samples is shown in Additional file 1: Figure S1.

Fig. 2

Transcriptional changes upon local infection by TYLCV in N. benthamiana. a Schematic representation of the experimental design. b Differentially expressed genes in the systemic infection are shown by volcano plot. The x-axis shows the log2 transformed gene expression fold change between infected and control samples. The y-axis indicates the negative log10 transformed adjusted p-values (FDR) of the differential expression test calculated by R package edgeR. The up-regulated and down-regulated genes are represented by red and blue dots, respectively. Pie chart shows the number of up−/down-regulated genes. c Validation of selected DEGs by qPCR. Values are the average of three biological replicates, relative to mock. NbACT was used as the normalizer. d Mapping of viral reads to the TYLCV genome. Three biological replicates (inner concentric circles) are represented; the upper side of each circle represents the virion (+) strand; the lower side of each circle represents the complementary (−) strand. ORFs are depicted in blue. The number of total reads are shown in Additional file 7: Table S2. e Validation of the expression of Rep and CP by qPCR. Expression values are relative to NbACT

A total of 247 and 1290 down- and up-regulated genes, respectively, were identified in these systemically infected samples (Fig. 2b; Additional file 8: Table S3). The RNA-seq results were validated by qPCR analysis of selected genes (Fig. 2c).

The clean reads were also mapped to the TYLCV genome, with an average of 3344 reads per million (RPM) (Fig. 2d; Additional file 7: Table S2). Also in this case, reads for both strands of the virus were detected throughout the viral genome with uneven and non-perfectly symmetrical distribution, and not restricted to the described ORFs; as observed in the locally infected samples, the number of reads was notably higher in the region of the genome containing the V2 and CP (late) genes (Fig. 2d; Additional file 7: Table S2). The accumulation of viral reads was confirmed by qPCR to detect expression of the viral genes encoding the Rep and the CP (Fig. 2e).

Distinct landscapes of transcriptional changes detectable in local and systemic infections by TYLCV in N. benthamiana

In both the locally and the systemically infected samples used, the virus is actively replicating, and therefore those putative transcriptional changes underlying successful viral multiplication must have been established. The vast difference in the number of differentially expressed genes (DEGs) between these infected samples raises the possibility that the larger proportion of infected cells in the agroinfected leaf patches might provide an increase in resolution, which may result from a lower dilution of infection-induced transcriptional changes due to a higher infected-to-uninfected cell ratio and/or negligible masking from potential non-cell-autonomous plant responses to the viral infection.

In order to gain further insight into the differences between locally and systemically infected samples, we set out to compare the subsets of induced or repressed genes in each case. Strikingly, as shown in Fig. 3a, both the up- and the down-regulated genes in systemic infections show only a partial overlap with those in local infections: only 56.7% of the repressed genes are also repressed in local infections, while, surprisingly, 6% are induced; among the induced genes in systemic infections, only 24% are also up-regulated in local infections, with a 23.8% down-regulated in these samples. Hierarchical clustering (Fig. 3b) shows that the systemically infected samples cluster closer to their control than to the locally infected samples. Taken together, these results clearly indicate that the differences detected between datasets go beyond a higher sensitivity in locally infected samples, which would explain the quantitative differences in DEGs, but not the apparent opposite behavior of some of them.

Fig. 3

Different transcriptional changes are detected upon local and systemic infections by TYLCV in N. benthamiana. a Venn diagram of DEGs in local and systemic infections. b The differential gene expressions of local and systemic infection samples are shown in heatmap with hierarchical clustering, which is drawn by R package pheatmap. The values in the gene expression matrix are centered and scaled in the row direction. The color bar indicates the Z-score. Loc: local infections; sys: systemic infections. EV: empty vector

With the aim of determining whether the distinct transcriptional landscapes of local and systemic infections may nevertheless result in similar functional outputs, we performed functional enrichment analysis using Gene Ontology and KEGG pathways annotations. As shown in Figs. 4 and 5 and Table 2, the overlap in over-represented GO categories (Biological Process Ontology) or KEGG pathways for systemic and local infections is only marginal. Common over-represented functional categories in both subsets of DEGs include trehalose biosynthetic process and defence response among the induced genes, and mitotic nuclear division and cellulose biosynthetic process among the repressed genes (Table 2).

Fig. 4

Functional enrichment analysis of differentially expressed genes in local and systemic TYLCV infections. Heatmap shows the Gene Ontology terms of biological processes that are significantly enriched (p < 0.01) in the up−/down-regulated genes in the local infection, the systemic infection, and their overlaps. Each square is colored according to the value of –log₁₀(p), where p is the p-value for the significance of GO term enrichment. The color bar indicates the –log₁₀(p). loc: local infection; sys: systemic infection; olp: overlapping between local and systemic infections. up: up-regulated; dn: down-regulated

Fig. 5

KEGG pathways of differentially expressed genes in local and systemic TYLCV infections. Heatmap shows the KEGG pathways that are significantly enriched (p < 0.05) in the up−/down-regulated genes in the local infection, the systemic infection, and their overlaps. Each square is colored according to the value of –log₁₀(p), where p is the p-value for the significance of KEGG pathway enrichment. The color bar indicates the –log₁₀(p). up: up-regulated; dn: down-regulated

Table 2 Significantly enriched GO terms (Biological Process Ontology) in the subsets of differentially expressed genes (p-value < 0.05)

Three KEGG pathways from Fig. 5 were selected to be graphically displayed, allowing for visualization and easy comparison of the transcriptional regulation of their components: MAPK signaling pathway (Additional file 3: Figure S3), Plant hormone signal transduction (Additional file 4: Figure S4), and Plant-pathogen interaction (Additional file 5: Figure S5). As observed in Additional file 3: Fig. S3 and Additional file 4: Figure S4, although not statistically significant in all cases, both local and systemic infections have an impact on these pathways, with the local infections having the strongest effect. As previously mentioned, the presence of the virus seems to activate plant defence responses (Additional file 3: Figures S3; Additional file 5: Figure S5). Notably, both types of infection trigger a detectable transcriptional repression of auxin signaling and a transcriptional activation of ethylene signaling, while only local infections resulted in a repression of the brassinosteroid signaling pathway (Additional file 4: Figure S4).

TYLCV and TbSCV modify the expression of a set of common genes upon systemic infection in N. benthamiana

In an attempt to identify potential central targets of the transcriptional geminiviral manipulation and/or effectors of the plant anti-geminiviral response, we decided to compare the transcriptional changes triggered by systemic infections by TYLCV and the geminivirus TbSCV, in combination with or without its associated satellite, in N. benthamiana [5]. Remarkably, a proportion of DEG were commonly affected by both TYLCV and TbSCV (12.6% of up-regulated genes by TYLCV infection, and 9.7% of down-regulated genes by TYLCV infection) (Fig. 6a); the proportion of induced, but not repressed, genes largely increased (to 28.7% of up-regulated genes by TYLCV infection) when TbSCV was inoculated in combination with its satellite (Fig. 6b), suggesting that some of the virulence functions provided by this ancillary molecule result in the activation of host genes and are already encoded in the TYLCV genome. Functional enrichment analysis of the genes commonly activated or repressed by TYLCV and TbSCV revealed the existence of a number of GO categories over-represented in these subsets (Table 3), including trehalose biosynthetic process and defence responses (in the common up-regulated gene set), and terpenoid biosynthetic process (in the common down-regulated gene set). Interestingly, a third geminiviral species, TYLCCN, has been proven to suppress terpenoid biosynthesis and release, and this effect in turn improves performance of its insect vector, the whitefly Bemisia tabaci [10]. Considering this negative regulation by three different geminivirus species, it is tempting to speculate that depletion of terpenoids is a requirement for geminiviruses to establish a successful infection in nature, perhaps at least partly through an indirect effect on favouring the insect vector-virus mutualism.

Fig. 6

Systemic infection by TYLCV and TbSCV in N. benthamiana trigger partially overlapping transcriptional changes. The transcriptional changes triggered by the systemic infection by TYLCV were compared to those triggered by the geminivirus TbSCV, without (a) or with (b) its associated satellite, in N. benthamiana; the TbSCV data are from Li et al., 2018

Table 3 GO enrichment of DEGs common to systemic infections by TYLCV and TbSCV in N. benthamiana. Over-represented GO categories (Biological Process Ontology) in the different subsets of DEGs represented in Fig. 6

Discussion

In this work, we describe and compare the genome-wide transcriptional changes detectable by RNA-seq occurring in N. benthamiana upon local or systemic infection by the geminivirus TYLCV. Our results show that, as expected and in agreement with previous works, infection by TYLCV causes a strong transcriptional reprogramming in the host; however, the detectable changes are more dramatic in our local infection system. It is possible that the local infection offers higher resolution owing to lower dilution of the infected cells; however, the opposite behavior of some DEGs in local and systemic samples suggests that the absence of transcriptional changes coming from non-infected cells and masking local changes also plays a role in the differences detected. Additionally, we cannot rule out that some differences between the two infected samples could be due to differences in the stage of the infection: although in both cases the virus is actively replicating in infected cells, in the leaf patch assay the infection is supposed to be more or less synchronized, while in the systemic infection a mixture of cells at different phases of the infection is most likely present in the samples. Another potentially relevant difference between local and systemic samples is the presence of Agrobacterium in the former; nevertheless, Agrobacterium is also present in the respective control. It is relevant to note that, if local and systemic samples show different results due to dilution and masking issues, these problems would potentially not only affect the study of transcriptional changes, but also the analysis of other molecular or physiological alterations resulting from the viral infection.

Functional enrichment analyses can help us gain insight into the biological processes underlying infection by TYLCV (Tables 2 and 3). In locally infected samples, over-represented GO categories included DNA recombination among the up-regulated genes, and cytokinesis among the down-regulated genes. These processes are expected to be connected to the viral manipulation of DNA replication and cell cycle, and might not be detectable in the systemically infected samples as a result of dilution. Photosynthesis also appears as transcriptionally down-regulated; photosynthetic shut-down seems to be a common outcome in viral infections ([8, 13,14,15,16], among others). A transcriptional negative regulation of BR signaling can also be detected in locally infected samples; these results are in agreement with a recent work by Seo et al. (2018), in which suppression of BR signaling was shown to underpin symptom development triggered by TYLCV in tomato.

Among those categories over-represented in both local and systemic infections, we can find cellulose biosynthetic process as transcriptionally repressed, suggesting that the viral infection might be affecting cell wall composition; changes in cell wall dynamics have been recently shown as triggered by Potato virus Y [17] and Rice tungro spherical virus [18], and in the latter case they have been proposed to correlate with virus-induced stunting. Intriguingly, both viruses negatively impact the cellulose biosynthetic machinery; whether impaired cellulose biosynthesis is a general plant response to the viral invasion is an idea that will require further investigation. Auxin signaling is also transcriptionally repressed upon the viral infection; these changes may also mediate or modulate the impact of the viral infection on plant development.

The identification of common subsets of up- and down-regulated genes in the systemic infections by TYLCV and TbSV [5] indicates that geminiviral manipulation of the cell and/or plant defence responses to geminiviral infection follow common transcriptional routes in different geminivirus/N. benthamiana interactions. Perhaps of particular interest is the finding that defence responses are activated in response to the viral infection; this observation reveals that these geminiviruses are being efficiently perceived as non-self by the plant, which in turn triggers a defence response. Although the activated plant defence reponses are not sufficient to fend off the virus, since the infection is established successfully, the fact that the plant is capable of detecting these pathogens in the first place paves the way for future engineering of the host, potentially boosting defences downstream of perception of the virus and hence tilting the balance in favour of the plant. Another pathway that emerges as a potential valuable target for engineering anti-geminiviral resistance is terpenoid biosynthesis, which is transcriptionally repressed by TYLCV and TbSV and has been proven to be suppressed by TYLCCN [10], raising the idea that its down-regulation might underpin a successful geminiviral infection.

Our results support the idea that the leaf patch assay entails great potential for the study of geminiviruses. Not only does this system result in a relatively synchronic infection and provides high resolution to detect virus-induced changes, but it also allows for the study of mutant viruses. All TYLCV null mutants for single genes, with the exception of those mutated in Rep, are capable of replicating their genome, but unable to infect the plant systemically; the use of local infections makes it possible to analyze the differences between the cellular changes triggered by wild-type and mutant viruses, therefore providing insight into the function of the viral genes in the context of the infection. An inherent limitation of this surrogate system, however, is the inevitable loss of cell type specificity: while TYLCV naturally infects phloem companion cells exclusively, in a leaf patch assay the virus is forced to replicate in mesophyll cells. Another obvious shortcoming is the impossibility of studying those mechanisms involved in cell-to-cell or long-distance transport, or the interactions between virus and vector.

All things considered, and while the analysis of systemic and local infections has provided and will continue to provide useful insight into the molecular events underlying the infection by geminiviruses, both approaches are imperfect for a number of reasons, as mentioned above. In order to considerably deepen our view, offering a substantial leap in our understanding of the molecular and physiological changes occurring during the plant-geminivirus interaction, isolating those infected cells, ideally based on the stage of the infection, will be crucial. Several approaches would enable the isolation of geminivirus infected cells: the use of transgenic plants harbouring a replicon-based system to label those cells sustaining active viral replication (like those described in [11]) could be combined with Fluorescence Activated Cell Sorting (FACS) or Laser Capture Microdissection (LCM), leading to the separation of infected from non-infected cells; high-throughput single-cell sequencing would also allow the unbiased identification and analysis of those cells containing the virus in the context of the infected plant. The increase in precision and resolution provided by the isolation of infected cells and their comparison to uninfected cells in the same plant will foreseeably result in an unprecedented view of the molecular landscaping triggered by the viral invasion.

Conclusions

Our results show that TYLCV induces a dramatic transcriptional reprogramming in N. benthamiana, the detection of which largely differs in local and systemic infections. Nevertheless, some responses, including a transcriptional repression of the auxin signaling pathway and a transcriptional activation of defence, can be commonly detected.

Comparison with the transcriptional changes induced by systemic infection by the geminivirus TbSV shows common subsets of up- and down-regulated genes similarly affected by both viral species, among which the suppression of terpenoid biosynthesis might be a general change triggered by geminiviruses. Taken together, our results not only provide insight into the transcriptional changes resulting from the infection by TYLCV in N. benthamiana, but also highlight the need to come up with an optimized system to gain a precise overview of the molecular and physiological changes caused in the host by the viral invasion.

Methods

Plant material and growth conditions

Wild-type N. benthamiana plants were grown in a controlled growth chamber in long day conditions (16 h light/8 h dark) at 25 °C.

Viral infections

The TYLCV infectious clone is described in [19, 20]; it contains a partial dimer of the TYLCV genome (AJ489258; [21]) in the pGWB501 vector [22]. Agrobacterium tumefaciens GV3101 strain was used for the delivery of TYLCV infectious clone and empty pGWB501 vector. Agrobacterium cells carrying these constructs were liquid-cultured in LB with appropriate antibiotics at 28 °C overnight. Bacterial cultures were centrifuged at 4000 g for 10 min and resuspended in the infiltration buffer (10 mM MgCl₂, 10 mM MES pH 5.6, 150 μM acetosyringone) to an OD₆₀₀ = 0.5. Bacterial suspensions were incubated in the buffer at room temperature and in the dark for 4 h before using them to infiltrate 4-week-old N. benthamiana for leave patch assays (local infections) and three-week-old N. benthamiana for systemic infection as described in [12].

RNA extraction

Total RNA was extracted from 8 mm leaf discs using the RNeasy plant mini kit (Qiagen) following the manufacturer’s instructions.

RNA sequencing

Transcriptome analyses were performed at the Genomic Core Facility, Shanghai Center for Plant Stress Biology, CAS. Three biological replicates were used. Total RNA (1 μg) from each sample was used for library preparation with NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England BioLabs, E7420L) following the manufacturer’s instructions. Prepared libraries were assessed for quality using NGS High-Sensitivity kit on a Fragment Analyzer (AATI) and for quantity using Qubit 2.0 fluorometer (Thermo Fisher Scientific). All libraries were sequenced in paired-end 125 bases protocol (PE125) on an Illumina HiSeq sequencer.

Quantitative RT-PCR

First-strand cDNA synthesis was performed with the iScriptTM cDNA Synthesis Kit (Bio-Rad #1708890) according to the manufacturer’s instructions. For qPCR reactions, the reaction mixture consisted of cDNA first-strand template, primers (500 nM each) and iTaqTM Universal SYBRR Green Supermix (Bio-Rad, #1725120). qPCR was performed in a BioRad CFX96 real-time system. Expression result was determined using the comparative Ct method (2^-ΔΔCt). Primers used are described in Additional file 9: Table S4. NbACT was used as the reference gene, using primers described in [23].

Preprocessing of RNA-Seq data

We cleaned the paired-end reads by Trimimomatic [24] (version 0.36). After trimming the adapter sequence, removing low quality bases and filtering short reads, clear read pairs were retained for further analysis.

Mapping and quantification of TYLCV reads

Cleaned reads were mapped to TYLCV DNA (GenBank: AJ489258.1) and its six ORFs by HISAT [25] (version 2.1.0) with default parameters. The RPM (Reads per Million) was used to quantify the expression level of each ORF and the whole viral genome. The read coverage of each base on the reference DNA was calculated by samtools [26] (version 1.5) with maximum coverage depth 8000 (−d 8000) and normalized to RPM. The expression level and read coverage were calculated for forward and reverse strand, respectively. The circular viral genome and read coverage of RNA-Seq data were visualized by CGView [27].

Reads mapping and quantification of N. benthamiana genes

The N. benthamiana draft genome sequence [28] (v1.0.1) was downloaded from the Sol Genomics Network (ftp://ftp.solgenomics.net/genomes/Nicotiana_benthamiana/assemblies/). Cleaned reads were mapped to the genome sequence by HISAT with default parameters. Number of reads that were mapped to each N. benthamiana gene was calculated with the htseq-count script in HTSeq [28].

Differential gene expression analysis

EdgeR [29] was used to identify genes that were differentially expressed between control and experiment samples. Genes with at least two-fold change in expression and had a FDR < 0.05 were considered differentially expressed genes (DEGs).

GO enrichment analysis

The Gene Ontology (GO) terms assigned to N. benthamiana genes were extracted from annotation file downloaded from the SGN (ftp://ftp.solgenomics.net/genomes/Nicotiana_benthamiana/annotation/Niben101/). GO enrichment analysis of DEGs was implemented by topGO [30] with custom gene-to-GOs mapping annotations.

KEGG pathway enrichment analysis

Since N. benthamiana is not supported by KEGG, we searched its ortholog genes in Arabidopsis by BLAST+ (version 2.5.0, evalue = 1e-5). The ortholog genes of the DEGs were used to perform KEGG enrichment analysis by clusterProfiler [31]. The enriched KEGG pathways were visualized with DEGs by Pathview [32].