Transcriptome analysis of Nicotiana benthamiana infected by Tobacco curly shoot virus

Background Tobacco curly shoot virus (TbCSV) is a monopartite begomovirus associated with betasatellite (Tobacco curly shoot betasatellite, TbCSB), which causes serious leaf curl disease on tomato and tobacco in China. It is interesting that TbCSV induced severe upward leaf curling in Nicotiana benthamiana, but in the presence of TbCSB, symptoms changed to be downward leaf curling. However, the mechanism of interactions between viral pathogenicity, host defense, viral-betasatellite interactions and virus-host interactions remains unclear. Methods In this study, RNA-seq was used to analyze differentially expressed genes (DEGs) in N. benthamiana plants infected by TbCSV (Y35A) and TbCSV together with TbCSB (Y35AB) respectively. Results Through mapping to N. benthamiana reference genome, 59,814 unigenes were identified. Transcriptome analysis revealed that a total of 4081 and 3196 DEGs were identified in Y35AB vs CK (control check) and Y35A vs CK, respectively. Both GO and KEGG analyses were conducted to classify the DEGs. Ten of the top 15 GO terms were enriched in both DEGs of Y35AB vs CK and Y35A vs CK, and these enriched GO terms mainly classified into three categories including biological process, cellular component and molecular function. KEGG pathway analysis indicated that 118 and 111 pathways were identified in Y35AB vs CK and Y35A vs CK, respectively, of which nine and six pathways were significantly enriched. Three major pathways in Y35AB vs CK involved in metabolic pathways, carbon metabolism and photosynthesis, while those in Y35A vs CK were related to Ribosome, Glyoxylate and dicarboxylate metabolism and DNA replication. We observed that 8 PR genes were significantly up-regulated and 44 LRR-RLK genes were significantly differentially expressed in Y35A treatment or in Y35AB treatment. In addition, 7 and 13 genes were identified to be significantly changed in biosynthesis and signal transduction pathway of brassinosteroid (BR) and jasmonic acid (JA) respectively. Conclusions These results presented here would be particularly useful to further elucidate the response of the host plant against virus infection. Electronic supplementary material The online version of this article (10.1186/s12985-018-1044-1) contains supplementary material, which is available to authorized users.


Background
Tobacco curly shoot virus (TbCSV), a monopartite begomovirus (genus Begomovirus, family Geminiviridae), was isolated from tobacco in Yunnan province, China in 2002 [1]. Some TbCSV isolates were identified associating with a betasatellite molecule (Tobacco curly shoot betasatellite, TbCSB) [2]. TbCSV is one of the most important pathogens causing leaf curl disease on tomato and tobacco, and is a severe constraint to crop yields. TbCSV consists of a circular single-stranded DNA genome that encodes 6 proteins, including AV1 (coat protein, CP) and AV2 in virion-sense strand whilst AC1 (replication-related, Rep), AC2 (transcriptional activator, TrAP), AC3 (replication enhancer, REn) and AC4 in complementary-sense strand [1]. TbCSB is a small circular single-stranded DNA molecule which encodes a sole protein βC1 [2]. Viral infection is a complicated procedure involving in the interaction between viruses and host plants. Understanding host responses to viral infection is advantageous in the development of effective strategies for virus control. In recent years, the interactions between viruses and host plants were clarified through analysis of transcriptomics [3][4][5]. RNA-seq is a recently developed approach to transcriptome profiling via deep-sequencing technologies, which provides a far more precise measurement of levels of transcripts and their isoforms than other methods [6]. Thus, this technique was widely applied in interpretation the interaction between virus infection and host plant. The multiple resistance mechanisms against cotton leaf curl disease (CLCuD) were revealed by the transcriptome analysis based on an RNA-seq in a naturally immune cotton species (Gossypium arboreum) caused by CLCuD [7]. The technology of transcriptomics and proteomics were employed to study the differential host gene expression during Mungbean yellow mosaic India virus (MYMIV) and Mungbean yellow mosaic virus (MYMV) infection under natural conditions [8]. The differential regulations of host genes revealed to be involved in cell cycle, cell-wall biogenesis, chloroplast, photosynthesis, hormone and sulphur assimilation pathways which may contribute to symptom development in soybean plants [8]. Recently, a research performed RNA-seq-based transcriptome sequencing of Jatropha curcas mosaic virus-infected and healthy leaf tissues of J. curcas, which provides a repertoire of molecular components after viral infection [9]. Similarly, several host genes were identified to be involved in different cellular and physiological processes during the infection of other viruses, including South African cassava mosaic virus (SACMV) [10], Tomato yellow leaf curl Sardinia virus (TYLCSV), Mungbean yellow mosaic India virus (MYMIV) [11,12], and Chilli leaf curl virus (ChiLCV) [13].
At present, it has been reported that the infectious clone of TbCSV induced severe upward leaf curling in N. benthamiana, but in the presence of TbCSB the symptom changed to a downward leaf curl [14]. Thus, TbCSV may represent an evolutionary intermediate between the truly monopartite begomoviruses and those that require the association of betasatellite [14]. Moreover, studies showed that AC2 and AC4 proteins of TbCSV mediate suppression of RNA silencing [15], and βC1 protein of TbCSB could bind single-and double-stranded DNA to suppress host RNA silencing activities [16]. However, it remains unclear that the molecular mechanism of viral pathogenicity, host defense and viral-host interactions. In this study, in order to get insights into the molecular mechanisms in response to TbCSV infection in N. benthamiana plants, the transcriptome analysis of N. benthamiana infected by TbCSV or TbCSV/TbCSB was performed, including differentially expressed genes, GO enrichment, KEGG enrichment, pathogenesis-related (PR) protein genes, LRR-RLK genes, brassinosteroid and jasmonic acid pathways. These data contribute to a better knowledge on the molecular mechanisms of TbCSV-and TbCSV/TbCSB-host interaction, and serve as a basis for devising new strategies to control TbCSV/TbCSB disease complex.

Virus sources and agroinoculation
The infectious clones of TbCSV isolate Y35 (pBinY35A-1.9) (Y35A) and its betasatellite (pBinY35β-2.0) (Y35B) were provided by Professor Xueping Zhou in the Biotechnology Institute of Zhejiang University. Both viral infectious clones were introduced into Agrobacterium tumefaciens strain EHA105. N. benthamiana (accession: Nb-1) plants were grown in a greenhouse under a 16 h light and 8 h dark cycle at 26°C. The Agrobacterium with infectious clones was infiltrated into N. benthamiana leaves at 4-6 leaf stage as previously described [2].

Materials for RNA sequencing
Leaf tissues of N. benthamiana plants were sampled at 20 dpi, and were immediately frozen and ground in liquid nitrogen and stored at − 80°C. After PCR detection, RNA-seq was performed for stored samples. Three treatments were setting as follows: plants infected with TbCSV alone (Y35A treatment); plants infected with TbCSV and TbCSB (Y35AB treatment); uninfected plants as control (CK). There are three biological replicates per treatment respectively.

RNA extraction and cDNA synthesis
Total RNA was extracted from N. benthamiana leaf tissues with RNAiso Plus (TAKARA Bio, Inc) as manufacturer's protocol. All RNA samples were first treated with gDNA Eraser (TAKARA Bio, Inc) and reverse-transcribed using a Prime Script RT reagent Kit (Perfect Real Time) according to manufacturer's instructions.

Library preparation for transcriptome sequencing
A total amount of 3 μg RNA per sample was used as input for RNA sample preparations. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purifed by using poly-T oligo-attached magnetic beads. After fragmentation, the cDNA was synthesized, and NEBNext Adaptors with hairpin loop structures were ligated to prepare for hybridization. PCR products were purifed (AMPure XP system), and library quality was assessed by Agilent Bioanalyzer 2100 system [18]. Finally, nine libraries were successfully constructed, and were sequenced using HiSeq™ 2500 equipment (Illumina, SanDiego, CA, USA) by 50 bp single-end (SE50) methods.

Quality control
For sequence quality control, raw data (raw reads) of fastq format were firstly processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. Meanwhile, Q20, Q30, GC-content and sequence duplication level of clean data were calculated. All the downstream analyses were based on the clean data with high quality [19].

Reads mapping to the reference genome
The draft sequence of the N. benthamiana reference genome have been downloaded from the SGN ftp site (ftp://ftp.solgenomics.net/genomes/Nicotiana_benthamiana/ assemblies) directly [20]. Index of the reference genome was built using Bowtie v2.2.3 (Broad Institute, Cambridge, MA, USA) [21] and single-end clean reads were aligned to the reference genome using TopHat v2.0.12 (Broad Institute, Cambridge, MA, USA) (mismatch = 2) [22]. The TopHat was selected as the mapping tool for that TopHat can generate a database of splice junctions based on the gene model annotation file and thus a better mapping result than other non-splice mapping tools [23].

Quantification of gene expression level
HTSeq v0.6.1 was used to count reads numbers mapped to each gene [24]. And then FPKM (Fragments per kilobase of transcript sequence per millions base pairs sequenced) of each gene was calculated based on the length of the gene and reads count mapped to this gene. FPKM, expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced, considers the effect of sequencing depth and gene length for the reads count at the same time, and is currently the most commonly used method for estimating gene expression levels [25].

Differential gene expression analysis
Differential expression analysis of Y35AB vs CK and Y35A vs CK was performed using the DESeq R package (1.18.0) (parameters: negative binomial distribution-based statistic, BH-FDR corrected p-value < 0.05) [26], which provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting p-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR).

GO term and KEGG pathway enrichment analysis
Gene Ontology (GO; http://www.geneontology.org) enrichment analysis of differentially expressed genes (DEGs) was implemented by the GOseq R package (corrected p-value < 0.05) [27], in which gene length bias was corrected. GO terms with corrected p-value less than 0.05 were considered significantly enriched by differential expressed genes. REVIGO was used for analysis of the enriched GO terms (http://revigo.irb.hr). KEGG (http://www.genome.jp/ kegg) is a database resource for understanding high-level functions and utilities of the biological system from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies. And KOBAS software was used to test the statistical enrichment of differential expression genes in KEGG pathways (corrected p-value < 0.05) [28].

Quantitative real-time PCR validation
To validate the RNA-seq data, quantitative real-time PCR (RT-qPCR) was conducted to examine the selected pattern of DEGs. For reverse transcription, 1 μg of total RNA was treated with gDNA Eraser (TAKARA Bio, Inc) and reverse-transcribed using a Prime Script RT reagent Kit (Perfect Real Time) and oligo (dT) as the primer according to the manufacturer's protocol. Primer sets were designed using AlleleID software (v6.0) (Additional file 1: Table S1). RT-qPCR was conducted by using NovoStart® SYBR qPCR Super Mix Plus (Novoprotein, Shanghai, China) on CFX 96 Real-Time System (Bio-Rad). Actin was used as a reference for calculating relative abundances using the 2 -△△Ct method [29]. All RT-qPCR experiments were performed in triplicate.

Plants with typical viral symptoms in Y35A and Y35AB treatments
At 20 dpi, Y35A inoculated plants showed slight upward leaf curling and Y35AB inoculated plants showed severe downward leaf curl (Fig. 1a-c). PCR detection showed that they are consistent with expectations ( Fig. 1d-e). Symptom observation and virus detection results showed that these plants were reliable for follow-up RNA-seq. Leaf tissues at same positions were collected from every treatment for RNA-seq.

Overview of transcriptome sequencing
To profile gene expression during virus infection, RNAseq libraries were constructed for the control-and virus-inoculated N. benthamiana plants. 5.76 GB clean bases data and 59,814 unigenes were obtained from nine samples.
There were 11834496-13118731, 12537134-15503750 and 11290015-12633203 clean reads in the three treatments of Y35AB, Y35A, and CK. The proportion of total reads that mapped to N. benthamiana genome of each sample was more than 96.09%, and that of uniquely mapped was above 72.69%. More than 95.58% of the clean reads had quality scores at the Q30 (an error probability for base calling of 0.1%) level [30]. Furthermore, the GC content ranges from 43.58 to 44.38%. The sequencing data are summarized in Table 1.
And 3196 transcripts (1391 up-regulated and 1805 downregulated genes) showed significant changes in Y35A vs CK (Fig. 2b, c and Additional file 3: Table S3). However, a total of 1232 DEGs in Y35AB vs CK overlap with those in Y35A vs CK, including 445 up-regulated and 786 down-regulated genes (Fig. 2c).

GO enrichment analysis
Gene Ontology, an internationally standardized gene function classification system, was used to classify the DEGs.
The results were enriched in Y35AB vs CK and Y35A vs CK, and 124 and 76 of them were significantly enriched, respectively. Further, the top 15 GO functional annotation terms were listed, and 10 of them were shared between Y35AB vs CK and Y35A vs CK, associated with biological process, metabolic process, cellular process, organic substance metabolic process, primary metabolic process, cellular metabolic process, nitrogen compound metabolic process, biosynthetic process, cellular nitrogen compound metabolic process and organic substance biosynthetic process (Fig. 3, Additional files 4 and 5: Tables S4 and S5). Furthermore, it was found that these enriched GO terms mainly classified into biological process.

KEGG pathway enrichment analysis
In order to further clarify molecular and biological functions of the genes, these DEGs were mapped to the KEGG database. KEGG pathway analysis indicated that 118 and 111 pathways were identified in Y35AB vs CK and Y35A vs CK, respectively, of which nine and six were significantly enriched with p-value< 0.05. The significantly enrichment analysis showed that DEGs of Y35AB vs CK involved in metabolic pathways, carbon metabolism, photosynthesis, carbon fixation in photosynthetic organisms, glyoxylate and dicarboxylate metabolism, porphyrin and chlorophyll metabolism, DNA replication, pentose phosphate pathway and nitrogen metabolism. While, DEGs involved in ribosome, glyoxylate and dicarboxylate metabolism, DNA replication, circadian rhythm-plant, photosynthesis-antenna proteins and nitrogen metabolism were significantly enriched in Y35A vs CK (Fig. 4, Additional files 6 and 7: Tables S6 and S7). Metabolic pathway was the major pathway which contains the largest number of DEGs in both Y35AB vs CK and Y35A vs CK, but DEGs involved in metabolic pathway were not significantly different in Y35A vs CK, it further indicated that the presence of TbCSB differentially changed metabolic pathways in N. benthamiana. The rich factor of porphyrin and chlorophyll metabolism pathway in Y35AB vs CK was near 1.00, the results showed that almost all genes in this pathway differentially expressed in TbCSV associated with betasatellite infected plants. Nevertheless, ribosome pathway might play an important role in response to TbCSV infection for the number of enriched DEGs in ribosome pathway wasnear to that enriched in metabolic pathways, and the rich factor was more than 0.75.
Screening of Y35AB vs CK and Y35A vs CK common highly expressed DEGs FPKM reflecting both the effect of sequencing depth and gene length for the read count, and is a commonly used method for estimating gene expression levels [31]. The fold change in gene expression was used to distinguish differentially expressed genes between samples because of the significance of digital gene expression profiles [32]. The intersection of highly expressed DGEs in Y35AB vs CK and Y35A vs CK was considered. Finally, 17 DEGs were identified (|log 2 (Fold change)| ≥ 1 and FPKM≥50) from significantly enriched KEGG pathways, these 17 genes were shared between Y35AB vs CK and Y35A vs CK all down-regulated, and involved in six pathways, including a b c photosynthesis, porphyrin and chlorophyll metabolism, metabolic pathways, carbon metabolism, carbon fixation in photosynthetic organisms, and glyoxylate and dicarboxylate metabolism. Additionally, six DEGs that not enriched in any pathway were screened (|log 2 (Fold change) | ≥ 2 and FPKM≥50), including four up-regulated and two down-regulated. In total, 23 candidate genes were screened as shown in Table 2.

Validation of the selected genes by RT-qPCR
To validate the RNA-seq data, 12 genes were selected randomly for RT-qPCR analysis from the Table 2, including genes GDCSa, GDCSb, PsaFa, PsaFb, PsaFc, PsaNa, PsaNb, glsF, CAX3, GRP, GRP3a and PU. The results indicated that the genes expression patterns by RT-qPCR were consistent with that by transcriptome sequencing (Fig. 5), the fold changes between RNA-seq and RT-qPCR were

Analysis of pathogenesis-related (PR) genes
An important feature of the plant defense in response to pathogen attack is the induction and accumulation of various PR proteins which are also a part of systematic acquired resistance (SAR). There were 29 PR genes identified by a local BLAST in RNA-seq dataset and online BLAST analysis based on reported PR family genes in N. tabacum. Among them, eight genes (27.6%) were significantly up-regulated, including PR-2 (Glucan endo-1,3-beta-glucosidase-like), PR-3 (Chitinase 8, glycoside hydrolase), PR-5 (Pathogenesis-related thaumatin superfamily protein), PR-11 (Chitinase-3-like protein 2), PR-17 (Plant basic secretory protein family protein), of them four genes were found in Y35AB vs CK, and five genes were observed in Y35A vs CK ( Table 3).

Analysis of LRR-RLKs genes
Leucine-rich repeat receptor-like protein kinases (LRR-RLKs) are the largest of receptor-like kinases in plants and play crucial roles in development and stress responses. In the present study, We combined RNA-seq dataset obtained LRR-RLK sequences with previously Fig. 5 Comparison of relative expression levels determined by RNA-seq and qPCR on 12 genes selected. All qPCR reactions were used three biological replica samples, repeated three times, and the data are presented as the mean ± SD. Relative expression levels are calculated from Ct values according to the 2 -ΔΔCt method. Actin was the reference gene for these qPCR experiments. Asterisks indicate statistically significant differences compared with CK, "*" indicate significant difference (0.01 ≤ p < 0.05), "**" indicate extremely significant difference (P < 0.01) NCBI reported sequences to perform analysis, results showed that a total of 71 unigenes were identified potential LRR-RLK genes in RNA-seq dataset, however, 44 significantly different unigenes were observed. Among them, 18 unigenes were found in Y35AB vs CK with 10 unigenes down-regulated and the other 8 unigenes up-regulated, and 32 unigenes expressed in Y35A vs CK, with 20 and 12 unigenes were down-regulated and up-regulated, respectively. Five unigenes were commonly regulated in both Y35AB vs CK and Y35A vs CK, including one up-regulated and four down-regulated unigenes (Table 4).

Analysis of brassinosteroid and jasmonic acid pathway
Plant hormones are not only essential for growth and development, but also play crucial roles in plant defense responses [33,34]. The transcriptome data showed that only one brassinosteroid (BR) synthesis-related gene, 3-epi-6-deoxocathasterone 23-monooxygenase (CYP90C1/ D1) was significantly decreased upon virus infection (Fig. 6a), it is likely result in a reduction of brassinosteroid. By contrast, expressions of jasmonic acid (JA) -ralated biosynthesis genes DOX1, HPL1, SAMT-X1 were decreased and PLA1/A2, loxF, KAT2, SAMT-X2 were significantly increased upon virus infection (Fig. 6b). We conjecture that the expression level of jasmonic acid might was increase. In addition, part of BR and JA signaling genes expression level were also altered, BRI1, BZR2, CYCD3-X2 were down-regulated and BZR1, CYCD3-X1 were up-regulated in BR signal transduction pathway, and JAR1, COI1, MYC2 were down-regulated and JAZ1, JAZ2 were up-regulated in JA signal transduction pathway (Fig. 6c). All synthesis-and signal transduction-related significantly difference expressed unigenes of BR and JA pathway shown in Table 5.

Discussion
TbCSV is a serious threat to many crops in China. In this study, RNA-seq -based transcriptomic characterization and comparative analysis of TbCSV-infected N. benthamiana with control plants was used to shed some light on understanding the molecular interactions of this pathosystem. 4081 and 3196 DEGs were screened from Y35AB vs CK and Y35A vs CK, respectively. Among the top 15 enriched physiological pathways, the porphyrin and chlorophyll metabolism, photosynthesis and pentose phosphate pathway were significantly affected in Y35AB vs CK. It suggested that, the disorder of these pathways might be responsible to the more serious leaf curling symptom on Y35AB-infected N. benthamiana plants. In addition, previous reports showed that βC1, an essential pathogenicity protein encoded by betasatellite, is required for symptom induction and RNA silencing suppression [2,16,[35][36][37]. Thus, βC1 protein might contribute to the specific changes of those genes or biological pathways in Y35AB-infected plants.
In our study, the DEGs expression with high levels that mainly concerned with photosynthesis-related genes including photosynthesis, porphyrin and chlorophyll metabolismand carbon fixation in photosynthetic organisms were down-regulated in N. benthamiana plants infected with TbCSV. Meanwhile, the KEGG pathway analysis showed significant enrichments of the photosynthesis, and the porphyrin and chlorophyll metabolism pathway in Y35AB-infected plants, with higher rich factor. It has been reported that photosynthesis-related genes could be affected upon virus infection such as Tobacco mosaic virus (TMV) [38], Cucumber mosaic virus (CMV) [39], Rice stripe virus (RSV) [40,41] and Alfalfa mosaic virus (AMV) [42]. And the chloroplast, one of the most dynamic organelles ofplant cell, plays an active part in the defense response and is crucial for inter organelle signaling [43]. Therefore, viruses need to suppress the chloroplast-mediated defense by employing pathogenicity factors, such as effector proteins, for successful infections [44,45]. Furthermore, a previous study showed that a Radish leaf curl betasatellite (RaLCB) damages the structural and functional integrity of chloroplasts, leading to inhibition of photosynthesis and symptom formation [46]. Some photosynthesis-related genes were down-regulated during both Y35A and Y35AB infections, indicating that these genes might participate in the interaction between TbCSV and N. benthamiana. An important feature of the plant defense response to pathogen attack is the induction and accumulation of various PR proteins, which are also a part of systematic acquired resistance (SAR) [47]. Previously, PR proteins were divided into five families based on their localization, isoelectric point, molecular mass, and biological activity [48]. PR proteins are currently categorized into 17 families. At present, besides PR-6 is predicted by automated computational analysis based on a genomic sequence (gene ID: LOC107768382), nine out of 17 families, including PR-1, PR-2 [49], PR-3, PR-4, PR-5 [50], PR-9 [51], PR-10 [52], PR-11 [53], PR-17 [54] families were reported that associated with the acquired resistance to pathogen infections. Among these PR proteins, PR-1 is generally considered as a marker gene of disease resistance [55]. Based on NCBI retrieval and local BLAST analysis, 29 PR genes in N. benthamiana were identified. In the present study, the expression levels of 8 PR genes were found to be significantly regulated in virus-infected N. benthamiana plants (Table 3). It was previously reported that in pathogens-infected wheat, the genes encoding peroxidase, PR-1, PR-2, PR-3, PR-4, and PR-5 were induced after 6-12 h and reached the highest levels at 36-48 h [56,57]. In our research, leaf tissues of N. benthamiana plants for RNA-seq assays were collected at 20 dpi, and the transcript levels of some PR genes (8 out of 29) were regulated significantly. For subsequent research, the expression of PR genes can be determined in a time course study. Based on our research, we speculated that some types of PR proteins might play roles in defense response to TbCSV infection.
In plants, many cellular signaling transductions are mediated by receptor-like kinases (RLKs). Leucine-rich repeat receptor-like protein kinases (LRR-RLKs) are the largest group of receptor-like kinases in plants and play crucial roles in development and stress responses [58]. LRR-RLKs regulate a wide range of plant growth and development, including meristem development [59,60], secondary growth [61], microsporogenesis and embryogenesis [62]. Besides, some LRR-RLKs mediated plant resistance against bacterial [63,64] or viral pathogens [65,66]. In this study, 71 LRR-RLKs genes were identified. Among them, 44 were differentially expressed in virus-infected N. benthamiana plants. The results showed that these LRR-RLKs genes might play important roles in N. benthamiana responses to TbCSV infection. By contrast, in the present of TbCSB, less number of LRR-RLKs genes were regulated (Table 4).
BR and JA are two important plant hormones that play crucial roles in plant growth and development [67] and also involved in plant defense responses to pathogen infections [68][69][70]. In the BR signaling pathway, BRASSINAZOLE RESISTANT (BZR) is a transcription factor, which orchestrating plant developmental and physiological processes by regulating BR-target gene expression [71,72]. BRASSINOSTEROID INSENSITIVE 1 (BRI1) is a brassinosteroid receptor, of which ubiquitously expresses plasma membrane-localized protein kinase [73], belonging to an LRR-RLKs subfamily [74]. For the JA signaling transduction, the perception of JA is mediated by the co-receptor CORONATINE INSENSITIVE 1 (COI1), which is an F-box protein cooperated with a family of the JASMONATE ZIM-domain (JAZ) transcriptional proteins. Upon the pathogens infection, bioactive JAs promote direct interaction between COI1 and JAZ proteins and then trigger the SCF-COI1 complex, resulting in poly-ubiquitination and degradation of the JAZ proteins. Ultimately, a variety of transcription factors could be activated for the expression of the JA-responsive genes [75]. Besides, MYC2 is a bHLH (basic helix-loop-helix) transcription factor and a positive regulator of JA responses [76]. A study has shown that rhizobacterium-mediated induction of JA reduces the symptoms of CMV infection in Arabidopsis thaliana ecotype Columbia plants [77]. In this study, BRI1, BZR1, BZR2, CYCD3, JAR1, COI1, JAZ and MYC2 genes were significantly regulated under the infection of TbCSV (Fig. 6). It suggested that BR and JA pathways played a role in response to TbCSV infection.

Conclusion
In this study, a genome-wide transcript profile was obtianed for N. benthamiana plant infected by TbCSV or plus TbCSB using RNA-seq approaches, and genes involved in plant defense system were found to be significantly regulated after TbCSV and TbCSV /TbCSB infection. The information provided in our study would be particularly useful for investigating the molecular mechanisms concerning begomovirus-host interaction and excavate resistance genes.