Symptoms and qRT-PCR experimental results after virus inoculation
After Xiangyan 11 was inoculated with HCRV and TSWV, the trends of the two viruses in pepper plants were inconsistent. The maximum accumulation of viruses was reached during the 1 dpi period. At 0–1 dpi, HCRV and TSWV both showed increasing trends.
During the period of 1–2 dpi, the accumulations of HCRV and TSWV gradually decreased; thus, 1 dpi was the time point for the maximum accumulation of the two viruses. During the time period of 3–7 dpi, the accumulation of HCRV and TSWV gradually increased. Finally, we chose 1 dpi samples and sent them to Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China) for high-throughput sequencing (Fig. 1).
Transcriptome sequencing results
To obtain the differences in gene expression of Xiangyan 11 inoculated with HCRV and TSWV viruses, we sequenced the leaves of Xiangyan 11 inoculated with the 2 viruses at 1 dpi on an Illumina HiSeq™ 2500 platform and constructed libraries with the sequencing data. In total, 12 cDNA libraries were prepared using poly(A) enrichment from 3 replicates of 4 treatments consisting of samples (NC1, NC2, NC3, CK1, CK2, CK3, T1, T2, T3 and H1, H2, H3). cDNA libraries underwent quality assessment and cleaning, and each sequencing library obtained more than 69 million clean reads (Supplemental Table S2). The Q20 and Q30 values of all libraries exceeded 98.46% and 94.76%, respectively, indicating that these data were suitable for further analysis (Table 1). The original data from the 12 cDNA libraries were saved to the SRA database with accession numbers SRR12726982- SRR12726994.
Table 1
Quality assessment of the cDNA libraries
Sample†
|
Raw Reads
|
Clean Reads
|
Clean Bases
|
Q20(%)
|
Q30(%)
|
GC(%)
|
NC1
|
87126070
|
85822070
|
12.77G
|
98.54
|
94.97
|
41.84
|
NC2
|
71022684
|
69909718
|
10.41G
|
98.46
|
94.76
|
41.91
|
NC3
|
74588286
|
73522478
|
10.95G
|
98.57
|
95.04
|
41.87
|
CK1
|
78260828
|
77222310
|
11.51 G
|
98.86
|
95.74
|
41.78
|
CK2
|
80078622
|
78869070
|
11.74 G
|
98.56
|
95.01
|
41.52
|
CK3
|
93153758
|
91816364
|
13.68 G
|
98.70
|
95.28
|
41.84
|
T1
|
79667856
|
78514498
|
11.68 G
|
98.52
|
94.90
|
41.86
|
T2
|
89897672
|
88601056
|
13.19 G
|
98.61
|
95.16
|
41.81
|
T3
|
90373984
|
89018518
|
13.24 G
|
98.53
|
94.96
|
41.83
|
H1
|
75646856
|
74693856
|
11.12 G
|
98.76
|
95.48
|
41.91
|
H2
|
72681602
|
71675354
|
10.67 G
|
98.55
|
95.01
|
41.81
|
H3
|
73708506
|
72669898
|
10.81 G
|
98.55
|
94.99
|
41.69
|
† NC indicates healthy plants, CK indicates plants inoculated with buffer, T indicates plants inoculated with TSWV, and H indicates plants inoculated with HCRV |
Clean reads were compared to the C. annuum reference genome (Capsicum Annuum.v.2.0, Accession number: AYRZ00000000) using Bowtie2 (2.2.8) and Tophat2 (2.1.1) software. The coverage of the comparison to the reference genome exceeded 92.30%, and the unique mapped reads exceeded 72.17%(Supplemental Table S3).
Cufflinks software was used to reconstruct transcripts based on Tophat2 comparison results and screen the positions of the assembled transcripts on the reference genome for new transcripts. The screening criteria were transcript length ≥ 200 bp and exon number ≥ 1. The number of known mRNAs obtained through screening exceeded 14,669, more than 34,434 new mRNA transcripts were obtained, and the total number of mRNA transcripts exceeded 49,181 (Supplemental Table S4).
We applied stringent filtering criteria to select highly significant differentially expressed genes (DEG) with at least two-fold differences between the two conditions by setting a cut-off q-value < 0.01 and log2-fold change > 1 (Fig. 2A). To better understand the diversity of DEG between each treatment, a Venn diagram was generated (Fig. 2B), and 36 DEG existed in all 3 comparison groups. The results showed that the expression levels of genes changed during HCRV and TSWV infection.
In this study, four comparison groups (CK-vs-H, CK-vs-T, T-vs-H and NC-vs-CK) were selected for functional annotation analysis. After obtaining the differentially expressed transcripts of each comparison group, GO (Gene Ontology) functional analysis and KEGG pathway (euKaryotic Orthologous Groups) analysis were performed on the DEG.
Compared with the control group (CK), the H treatment group recorded a total of 2617 DEG belonging to 40 functional groups, and the T treatment group recorded a total of 1587 DEG belonging to 37 functional groups, including biological processes, cell components and molecular functions (Supplemental Figure S1).
To further clarify the response mechanisms of DEG to virus infection of plants, the KEGG database was used to classify and characterize transcriptome DEG into corresponding pathways. Compared with the control group, the H comparison group (CK-vs-H) had 658 enriched DEG in 115 pathways. Among these pathways, the 13 pathways with the most significant enrichment (p-value < 0.05) were ribosome, protein processing in the endoplasmic reticulum, carbon fixation in photosynthetic organisms, photosynthesis-antenna protein, stilbenoid, diarylheptanoid and gingerol biosynthesis, pyruvate metabolism, biosynthesis of amino acids, flavonoid biosynthesis, alanine, aspartate and glutamate metabolism, DNA replication, nitrogen metabolism, microbial metabolism in diverse environments, and arginine biosynthesis. In addition, metabolic pathways related to plant immunity, plant-pathogen interactions, plant hormone signal transduction and other pathways were involved in plant immunity DEG (Fig. 3A). Compared with the control group, the T treatment group (CK-vs-T) exhibited 315 enriched DEG in 103 pathways. Among these pathways, the 14 pathways with the most significant enrichment (p-value < 0.05) were DNA replication, phenylpropanoid biosynthesis, stilbenoid, diarylheptanoid and gingerol biosynthesis, biosynthesis of secondary metabolites, flavonoid biosynthesis, spliceosome, α-linolenic acid metabolism, pyruvate metabolism, metabolic pathways, plant circadian rhythm, amino sugar and nucleotide sugar metabolism, riboflavin metabolism, glutathione metabolism, and protein processing in the endoplasmic reticulum (Fig. 3B). In the T-vs-H comparison group, 251 DEG were enriched in 99 pathways. Among these pathways, the five pathways with the most significant enrichment (p-value < 0.05) were the photosynthesis-antenna protein pathway, ribosome, glutathione metabolism, carotenoid biosynthesis, and monoterpenoid biosynthesis (Fig. 3C). In the NC-vs-CK comparison group, 125 DEG were enriched in 81 pathways. The 10 pathways that were significantly enriched (p-value < 0.05) were photosynthesis-antenna protein pathways, diterpenoid biosynthesis, pyruvate metabolism, amino sugar and nucleotide sugar metabolism, biosynthesis of unsaturated fatty acids, phenylalanine, tyrosine and tryptophan biosynthesis, tryptophan metabolism, α-linolenic acid metabolism, cutin, suberine and wax biosynthesis, and biosynthesis of secondary metabolites (Fig. 3D).
Identification and analysis of plant responses to DEG
To further analyse the immune responses of the two viruses in infected plants, we focused on the plant-pathogen interaction pathways, plant hormone signal transduction and phenylpropanoid synthesis pathways enriched by DEG (Fig. 4, Supplemental Table S5).
When a plant is infected by a pathogen, it triggers its own defence response and usually uses two defence mechanisms to resist invasion of the pathogen. In total, 18 genes were differentially expressed in the comparison group of H and control (CK-vs-H), 11 genes were differentially expressed in the comparison group of T and control (CK-vs-T), and 2 genes were differentially expressed in the comparison group of T-vs-H. Therefore, we found that the defensive response of HCRV to Xiangyan 11 was stronger than that of TSWV.
In the XY11-H sample, we found that most DEG were significantly upregulated (8 upregulated genes, 2 downregulated genes), including calcium-dependent protein kinase (CDPK), heat shock protein (HSP90), and LRR receptor protein kinase (FLS2), with the highest expression at 2.05. Heat shock proteins play a role in the defence process of ETI, and these genes cause a hypersensitive response (HR) in the process of plant resistance to pathogen infection (Fig. 4A).
In the plant hormone signal transduction pathway, it was found that many hormones were upregulated after infection with TSWV and HCRV. HCRV cannot successfully infect Xiangyan 11, but TSWV can successfully infect Xiangyan 11. After HCRV inoculation with Xiangyan 11, the expression levels of TIR1 and SAUR36 genes in the auxin synthesis pathway were significantly downregulated. The upregulated expression of SAPK1 in the abscisic acid (ABA) pathway may affect the synthesis of abscisic acid. In the synthesis of ethylene (ETH), EBF1 downregulates its expression, which relieves the inhibition of the ubiquitinated proteolytic pathway. In the jasmonic acid (JA) synthesis pathway, the two TIFY10A genes in the zinc jasmonic acid domain (JAZ) and the transcription factor MYC2 (bHLH14) are upregulated, and the upregulated expression of JAZ promotes the process of ubiquitinated proteolysis. The upregulated expression of JAZ and MYC2 is conducive to the synthesis of JA and the defence response of plants. The key genes histidine kinase 2 (AHK2) and B-ARR (ARR12) in the cytokinine synthesis pathway (CTK) are both upregulated, and the expression of A-ARR (ARR9) is downregulated. Since the A-ARR gene inhibits the expression of the B-ARR gene, the upregulated expression of the abovementioned genes ultimately increases the accumulation of CTK (Fig. 4B).
The secondary metabolism of plants plays an important role in the defence response of plants, and the synthetic pathway of phenylpropanoids is a common mechanism. The final product of phenylpropanoid biosynthesis is lignin. The accumulation of lignin on the cell wall forms a physical barrier to pathogenic bacteria. In the phenylpropanoid biosynthesis pathway, the expression levels of caffeic acid methyltransferase (COMT), anti-cinnamic acid 4-monooxygenase (CYP73A), peroxidase (PER) and oxalate O-hydroxycinnamyl transferase (HST) were upregulated after HCRV and TSWV infection. When HCRV infected peppers, acetyl-CoA-benzyl alcohol acetyltransferase (BEAT), cinnamyl alcohol dehydrogenase (CAD), peroxidase (PER), oxalate O-hydroxycinnamyl transferase (HST), and CYP73A were all upregulated, increasing the accumulation of lignin, which is beneficial to the defence response of plants (Fig. 4C).
Changes in transcription factor expression in the plant defence response
Transcription regulation is an important mechanism of eukaryotic gene expression and plays an important role in plant defence responses. Transcription factor family (TF) analysis of DEG in 12 cDNA libraries showed that there were 2350 DEG distributed in 52 transcription factor families, of which C3H (674), bHLH (142), MYB (130), C2H2 (111), ERF (109), MYB-related (103), NAC (77), and WRKY (72), among others, are related to plant defence response (Fig. 5).
Transcriptome data showed that compared with the control (CK), a total of 32 TFs were expressed after HCRV treatment, including 17 upregulated and 15 downregulated TFs. These TFs are divided into 16 different families, mainly C3H, MYB, ARR-B, ERF, LBD, HD-ZIP, bHLH and MYB-related (Fig. 6A). Among them, the upregulated TF families mainly included MYB (12.5%), C3H (9.4%), ARR-B (6.3%), bHLH (6.3%) and MYB-related (6.3%), and the downregulated TFs mainly included ERF (9.4%), C3H (6.3%), ARR-B (6.3%) and LBD (6.3%).
Compared with the control (CK), 33 TFs were differentially expressed after TSWV treatment, including 20 upregulated and 13 downregulated TFs. These TFs are divided into 21 different families, mainly ERF, HD-ZIP, MYB-related, bHLH, CO-like, WRKY, ARR-B, and AP2 (Fig. 6B). Among them, the upregulated TF families mainly included MYB-related (6.1%), bHLH (6.1%), CO-like (6.1%), and HD-ZIP (6.1%), and the downregulated TF families mainly included ERF (9.1%), ARR-B (6.1%), and LBD (6.1%). Many of these TF are involved in the response to plant defence and signal pathways activated by JA. Therefore, we focused on the transcription factors WRKY, MYB, AP2/ERF and other reactions to plant defence.
We analysed the TFs in the CK-vs-H and CK-vs-T comparison groups. The results showed that TFs such as C3H, MYB, and ARR-B in the CK-vs-H comparison group were upregulated, but WRKY, ERF and LBD TFs were downregulated (Fig. 6A); in the CK-vs-T comparison group, MYB, MYB_related, bHLH, bZIP and other TFs were upregulated, and ERF, WRKY and other TFs were downregulated (Fig. 6B).
qRT-PCR experiment results
To further confirm the gene expression pattern of the differentially expressed transcripts identified by RNA-Seq, we selected 10 annotated unigenes for qRT-PCR experiments (Fig. 7). These candidate genes are closely related to the transcription factor family and hormone and phenylpropanoid synthesis. The results showed that the 10 gene expression profiles obtained by qRT-PCR had high similarities with RNA-Seq data. These findings confirm the reliability of the differential expression analysis data obtained by RNA-Seq and reflect the actual transcriptomic changes in the defence response of peppers after virus infection.