The tobacco varieties used in this study were CBH (RKN-susceptible) and G28 (RKN-resistant) (Zhu et al., 1989; Li et al., 2010), which were obtained from the breeding laboratory of Henan Agricultural University (Henan, China). PVY M^SN^R was obtained from the Yunnan Academy of Tobacco Agricultural Sciences (Yunnan, China). Tobacco seedlings grown to five or six true leaves were inoculated with PVY M^SN^R solution (comprising 40 ml ddH2O and 1 g leaves of preserved PVY M^SN^R diseased plants), and the leaf surface was lightly rubbed with arenaceous quartz to cause minor wounds. After infection, plants were placed in growth incubators with a photoperiod of 12 h light/12 h dark cycle, light intensity of 1500 lux, with temperatures of 26°C in the light and 22°C in the dark. Tobacco leaves before infection were set as the control group, and tobacco leaves at 1, 3, 5, and 7 days after infection were set as treatment groups; each treatment was repeated three times, and each replicate was prepared by mixing leaves from five plants in each treatment group.

To measure H₂O₂, after infection for 7 d with PVY M^SN^R, leaves were incubated in 1 mg ml^–1 diaminobenzidine (DAB) solution (pH 3.8) at room temperature for 8 h. Stained seedlings were then transferred to 70% (v/v) ethanol to remove the chlorophyll. DAB staining was repeated three times.

RNA isolation was performed using the mirVana miRNA Isolation Kit (Ambion, Inc., Austin, TX, USA). To determine the degradation, contamination, concentration, and integrity of RNA, analyses using 1% agarose gels, a Qubit® RNA Assay Kit with a Qubit®2.0 Fluorometer (Life Technologies, CA, USA), and a Bioanalyzer 2100 system (Agilent Technologies, CA, USA) were conducted. Subsequent analysis was performed using samples with an RNA Integrity Number ≥ 7. Library construction was performed using a TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, USA). Thirty libraries were sequenced by Shanghai OE Biotech on the Illumina NovaSeq platform (Illumina, San Diego, CA, USA) to generate 150-bp paired-end reads. All sequence data were deposited in the NCBI database (BioProject accession: PRJNA984343).

Trimmomatic version 0.36 was used to process raw reads (raw data) (Bolger et al., 2014). Low-quality reads and poly-N reads were filtered out, and high-quality, clean reads were finally obtained. Clean reads were mapped to the reference genome using Hisat2 version 2.2.1.0 (Kim et al., 2015).

Differences between G28 and CBH transcriptome data were assessed by principal component analysis (PCA). Cufflinks version 2.2.1 was used to calculate fragments per kilobase per million reads (FPKM) values (Trapnell et al., 2010) and htseq-count to obtain read counts for all genes (Anders et al., 2015). Expression trend analysis was performed using OmicShare Tools (https://www.omicshare.com/tools/Home/Soft/trend); FPKM data for genes in each sample were imported, and the default settings were used. The R package, DESeq version 1.18.0, was used to identify DEGs. Significantly differential expression was defined as P < 0.05 and fold-change > 2 or < 0.5 (Anders and Huber, 2010). The Venn and column diagrams were drawn using TBtools with default settings (Chen et al., 2020). R packages based on the hypergeometric distribution (GO and KEGG) were used for gene function analysis (Kanehisa et al., 2008). The top 20 functional pathways of each group were selected for Heatmaps plotting. The heatmaps of functional analysis were based on Log₁₀ ^{(p-value of each pathway)} through OmicShare Tools with the default settings.

The changes in expression of six randomly selected genes were confirmed by conducting qRT-PCR using an ABI 7500 fast Real-Time PCR System (Applied Biosystems, Waltham, MA, USA). Briefly, the total RNAs were isolated as described previously. A TranScript All-in-One First-Strand cDNA Synthesis SuperMix for qPCR Kit (Transgen Biotech, Beijing, China) was employed to synthesize the cDNA. qRT-PCR was performed following the protocol of the TransStart Top Green qPCR SuperMix kit (Transgen Biotech, Beijing, China). The Primer-BLAST online NCBI tool was used to design primers specific for each DEG (Table S1). The L25 gene was used as the internal reference to normalize the relative expression levels. The gene expression level of C0 (before infection of CBH) and G0 (before infection of G28) were set to one, and the 2^−ΔΔCt method was used to calculate the relative expression level of each assessed gene (Livak and Schmittgen, 2001). Sequencing data were verified by comparing the expression trends of each gene as determined by qRT-PCR and FPKM values. Three biological technical replicates were performed to validate the results of qRT-PCR.

Leaves of the two tobacco varieties were compared after infection for 7 d with PVY M^SN^R. As shown in Figures 1A–D, the RKN-resistant variety, G28, developed necrosis spots caused by HR at 7 d, while the RKN-susceptible variety, CBH, did not. Qualitative analysis of hydrogen peroxide in the two varieties after infection was conducted by DAB staining (Figures 1E–L). Hydrogen peroxide content did not change significantly in CBH leaves after infection; however, hydrogen peroxide staining was clearly observed in G28 leaves one day after infection and increased and became more diffuse in the following days.

RNA sequencing data statistics are presented in Table S2. The overall gene expression levels of the three biological replicates of each sample were similar, and the identified genes showed differences in expression levels on different days after infection in the two varieties (Figure 2A). The repeatability and differentiation of each sample were assessed by PCA analysis. The three replicates of each sample clustered together, and G28 and CBH samples were easy to distinguish on the PC1 axis (Figure 2B).

Expression trend analysis showed that three and seven significantly enriched profiles (p < 0.05) were identified in G28 and CBH samples, respectively. Two model profiles (profile 0, downregulated profile; profile 19, upregulated profile) were identified as common to both varieties (Figure 3). The potential functions of genes in profiles 0 and 19 in the two varieties are shown in Figure S1.

Numbers of DEGs 1, 3, 5, and 7 d after infection, each compared with 0 d, were calculated for both plant varieties. More up-regulated DEGs were found in CBH at each stage, while more genes were suppressed in G28. The number of DEGs between the two varieties at different time points is also presented (Figure 4A).

After infection, the number of common upregulated and downregulated DEGs was relatively higher in G28 (1347 and 3092 DEGs, respectively) than in CBH (1197 and 894 DEGs, respectively) at 1, 3, 5, and 7 d compared with normal conditions (Figure 4B).

The results of KEGG pathway enrichment analysis for both varieties after infection for 1, 3, 5, and 7 d compared with uninfected controls are shown in Figure 5. The ‘Photosynthesis - antenna protein’, ‘Metabolism of xenobiotics by cytochrome P450’, and ‘Drug metabolism - cytochrome P450’ pathways were common to both G28 and CBH varieties after inoculation. Further, ‘Plant hormone signal transduction’ and ‘Phenylpropanoid biosynthesis’ were specifically enriched in CBH after inoculation, while ‘Glyoxylate and dicarboxylate metabolism’, ‘Glutathione metabolism’, ‘Biosynthesis of amino acids’, and ‘Carbon metabolism’ were specifically enriched in G28 after infection. Significant functional differences among DEGs between G28 and CBH at each stage after infection are shown in Figure 5C. Photosynthesis-related metabolic pathways (‘Porphyrin and chlorophyll metabolism’, ‘Carbon fixation in photosynthetic organisms’, ‘Photosynthesis’, and ‘Photosynthesis - antenna proteins’), ‘Plant hormone signal transduction’, ‘Starch and sucrose metabolism’, and ‘Phenylpropanoid biosynthesis’ were the main metabolic pathways that differed between the two varieties.

GO enrichment analysis was also performed to identify the putative function of the DEGs after infection (Figures S2–S4).

DEGs involved in abscisic acid (ABA) signal transduction were identified. Further, four and three PYL genes were highly expressed in CBH and G28 leaves, respectively. Notably, all identified ABF, ABI, and PP2C-like genes were expressed at higher levels in G28, and the expression levels of SnRK genes showed opposite trends in the two varieties (Figure 6A). DEGs involved in the brassinosteroid (BR) signal transduction pathway were obtained, of which CYCD3, BKI1, and one gene were highly expressed in CBH; two BSK genes (BSK8-1 and BSK8-2) were highly expressed in G28; and BZR2-1 was highly expressed in CBH, while BZR2-2 was highly expressed in G28 (Figure 6B). DEGs involved in the cytokinin (CK) signal transduction pathways were obtained; almost all ARR genes were highly expressed in CBH (Figure 6C). All identified DEGs in the ethylene, gibberellin, and salicylic acid (SA) signal transduction pathways were highly expressed in G28 (Figures 6D–F). Moreover, most indole-3-acetic acid (IAA) signal transduction pathway genes were highly expressed in CBH (Figures 6G, H).

In general, most SUS and CWINV genes were expressed at higher levels in G28 than in CBH samples (Figure 7), while levels of SUS1, SUS4, and SUS5 showed the opposite trend in the two varieties. SUS2, SUS3, CWINV4-1, and CWINV4-3 expressions gradually decreased after infection in both varieties. CWINV4-2 gene expression also presented a reduced trend in CBH but began to increase at 3 d after infection in G28, reaching a peak at 5 d after infection. The BGL, PGI, and PGMP genes were all highly expressed in CBH relative to G28 leaves. Most BGL, PGI, and PGMP genes showed an increased trend in CBH but decreased continuously in G28. Expression levels of four SPS genes presented a reduced trend in G28 leaves but gradually increased and reached peak expression levels at different stages after infection in CBH. Two HXK genes (HXK2-1 and HXK2-2) were downregulated after infection in CBH, with minimum levels at 1 d. In contrast, both HXK genes were upregulated after infection in G28, with the maximum levels at 7 and 1 d, respectively. One FPK gene was downregulated after infection in CBH, with minimum levels at 7 d; however, its expression was upregulated after infection in G28, with the highest levels at 5 d. Four AGPS genes were identified and increased after infection in CBH with maximum levels at 1 d (AGPS1-1, 1-2, and 1-3) and 5 d (AGPS1-4). AGPS gene levels decreased after infection in G28, with the lowest levels at 3 d (AGPS1-1 and 1-2) and 7 d (AGPS1-3 and 1-4). Levels of all six SS genes were similar to those of AGPS in the two varieties. Two SBE genes (SBE1 and SBE2) decreased after infection in the two varieties but were expressed at higher levels in CBH than in G28 samples. SBE3 was also downregulated after infection in CBH but showed an increased trend in G28, with maximum levels at 5 d.

Fifty-six DEGs involved in phenylpropanoid biosynthesis were identified (Figure 8). All of the six PAL and three 4CL genes were expressed at higher levels in G28 than in CBH strains, particularly before infection. PAL gene expression levels gradually decreased after infection in CBH, while in G28, PAL1, PAL2, PALA1, and PALA2 presented maximum expression levels before infection, with PAL3 and PAL6 reaching peak levels at 5 d. Levels of the 4CL genes gradually decreased after infection in both CBH and G28 plants. In CBH, four CAD genes were induced suddenly after infection and presented maximum expression levels at 1 d, followed by a decline. CAD genes gradually decreased after infection in G28, except for CAD1-2, which showed the highest expression level at 5 d.

The levels of DEGs associated with photosynthesis-related metabolic pathways are shown in Tables S3–S6. Ninety-one DEGs related to photosynthesis were identified, most of which were induced after infection in CBH, and had maximum expression levels at 1 d, followed by 5 d. These genes were expressed at lower levels in G28 after infection, with minimum levels from 3 to 7 d after infection. Sixty DEGs related to antenna proteins were identified. Most genes were reduced at 1 d, then increased at 3 d after infection and reached peak levels at 5 or 7 d in CBH. In G28, levels of these genes decreased significantly after infection and remained low until 7 d. Seventy-three DEGs involved in carbon fixation in photosynthetic organisms were obtained, of which most levels were reduced after infection in CBH and decreased significantly after infection in G28. Fifty-one DEGs related to porphyrin and chlorophyll metabolism were identified; the majority were induced after infection and reached maximum levels at 1 d, then declined in CBH. In G28, expression of these genes was significantly reduced after infection and remained low for 7 days.

Six DEGs were randomly selected for validation of sequencing data by qRT-PCR. The qRT-PCR assay results exhibited general and high consistency with the sequencing data, indicating that our RNA-seq analysis findings are reliable (Figure S5).