The seeds of petunia (Petunia hybrida) cultivar ‘Primetime Blue’ were purchased from Goldsmith Seeds Inc. (Gilroy, CA, United States). They were planted in a 72-well plastic tray filled with the soil mix of peat moss, perlite, and vermiculite at 2:1:1 by volume. The tray was placed into a growth chamber at 25/22°C temperature (day/night) with a photoperiod of 16/8 h (light/dark) and a 70% relative humidity. Young seedlings were transplanted to small plastic pots containing the same soil mix for inoculation with wild TRV strain and Agrobacterium culture bearing TRV-GFP-PhPDS or pCAMBIA2300 constructs. Wild TRV strain-inoculated petunia leaves at three stages were used for Illumina RNA-Seq. The leaves infected with various TRV constructs were used for VIGS assay. The leaves infiltrated with pCAMBIA2300 constructs were used for transient overexpression assay. To examine the abundances of transcripts, the inoculated or uppermost systemically infected leaves were collected for RNA extraction.

Tobacco rattle virus strain PPK20 was preserved in tobacco (N. tabacum) cultivar ‘NC82’ plantlets. For infectious sap preparation, the young leaves infected with TRV (PPK20) were harvested and homogenized with a small amount of kieselguhr in 100 mM phosphate buffer at pH 7.0 (1:6, w/v) using a mortar and pestle. The viral preparations were mixed with 10% of 600-mesh carborundum powder by weight, and then mechanically inoculated onto the fully expanded healthy young leaves (Tian and Valkonen, 2013). To ensure sufficient viral propagation, the rub-inoculation process was performed once more within half an hour. The inoculum was flushed with sterile deionized water. The symptom occurrence was observed during the post-inoculation periods. Relative accumulation levels of TRV (RNA1 and RNA2) were determined by real-time quantitative PCR (RT-qPCR) analysis.

Petunia leaves infected with TRV (PPK20) were collected to extract total RNA using TRIzol reagent (Invitrogen, Carlsbad, CA, United States). Three independent biological replicates were used for each stage. RNA samples were further purified to eliminate contaminant DNA using RNase-free DNase I (Promega, Madison, WI, United States). RNA integrity was examined by observing the major RNA species (18S and 28S rRNA) through 1.5% agarose gel electrophoresis. RNA yield was quantified using a NanoDrop ND-2000c Spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States). Illumina sequencing and generation of transcriptome data were conducted at Gene Denovo Bio-Technology Co., LTD. (Guangzhou, China). In brief, the general experimental procedure includes mRNA enrichment and fragmentation, cDNA synthesis, adaptor ligation, and ligated products sequencing as previously described (Sun et al., 2019).

High-quality clean reads were produced by removing adaptor, low-quality, and unknown nucleotides of raw reads through a Perl script tool. The resulting clean reads were assembled via Trinity program (v2.1.1) (Grabherr et al., 2011). The redundancy of assembled products was removed to obtain unique fragments using iAssembler software (Zheng et al., 2011). The fragments over 200 bp were enriched using TIGR Gene Indices clustering tools (TGICL, v2.1) (Pertea et al., 2003), and regarded as unigenes. The unigenes were mapped to the reference genome¹ of P. axillaris, which is one of P. hybrida’s parents. Novel transcripts were determined using Cuffcompare software through a threshold for length ≥ 200 bp and exon number ≥ 2 (Trapnell et al., 2012), and annotated against the public protein databases including NR, Swiss-Prot, COG, and KEGG.

Read counts of unigenes were calculated and standardized to the fragments per kilobase of transcript per million mapped reads (FPKM), which were used to quantify the expression levels of unigenes through Cufflinks program (v2.1.1) (Trapnell et al., 2010). Differentially expressed genes (DEGs) were identified using edgeR package (v3.12) (Chen et al., 2014). Statistical significance of DEGs was evaluated with a cutoff of fold change (FC) ≥ 2.0 and false discovery rate (FDR) ≤ 0.05. Blast2GO tool (v2.3.5) and WEGO program² were used to gain knowledge about Gene Ontology (GO) classification of DEGs (Conesa et al., 2005; Ye et al., 2006). To assign putative metabolic processes to DEGs, the pathway enrichment analysis was carried out using a BLAST search against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database³ (Wixon and Kell, 2000). P-values were calculated using a previously reported method and subjected to FDR correction (Lu et al., 2012). Significantly enriched pathways were characterized using FDR ≤ 0.05 as a cutoff.

Petunia leaves inoculated with TRV (PPK20) at different time points were harvested. For ethylene (ET) detection, the leaves were placed in a 50-ml plastic tube at 25°C for 4 h, with an airproof rubber plug embedded into the cap. The parafilm was used to seal the connection gaps. A 3-ml gas sample was withdrawn from the tube using a gas-tight syringe with a needle, and applied to a gas chromatograph (GC-8A; Shimadzu, Kyoto, Japan) for measurement of ET release. The examination of other hormones, including abscisic acid (ABA), jasmonic acid (JA), and salicylic acid (SA), was performed as previously described (Mornya and Cheng, 2011) with some minor modification. Briefly, petunia samples were extracted in 80% methanol and 1 mM butylated hydroxytoluence on a rotary shaker at 4°C for 1 h. After a complete homogenization, the mixture was centrifuged at 12,000 rpm for 20 min at 4°C. The supernatant was collected and purified through a 0.45 μm filter. It was concentrated to near dryness in N₂ gas and then dissolved in 50% methanol. The resulting solution was analyzed by a high performance liquid chromatography-electrospray ionization tandem mass spectrometry (HPLC-ESI-MS/MS). All chemical substances purchased from Sigma-Aldrich (St Louis, MO, United States) were used as the reference standards.

Total RNA from petunia leaves was extracted as mentioned above. After purification and quantification, about 2–5 μg of RNA samples were used to synthesize first-strand cDNA with SuperScript III reverse transcriptase (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s recommendations. The resulting cDNA products were fivefold diluted by ddH₂O and restored at –20°C until further use. The specific primer pairs were designed based on the known cDNA sequences of genes tested via Primer Premier5⁴ (Supplementary Table S1). 26S ribosomal RNA served as a reference gene. RT-qPCR assay was performed using the SYBR Green PCR Master Mix (2×) (Applied Biosystems, Foster City, CA, United Stataes) in a LightCycler480 Real-Time PCR System (Roche Diagnostic, Basel, Switzerland). The experiment was repeated three times using independent RNA samples. Relative expression levels of genes were analyzed using a previously described protocol (Livak and Schmittgen, 2001).

The TRV-GFP construct was generated as previously described (Tian et al., 2014). A 138-bp fragment of PhPDS cDNA was amplified using the forward and reverse primers with a SacI and an XhoI restriction site, respectively. It was ligated into the corresponding sites of TRV-GFP plasmid to form the TRV-GFP-PhPDS construct. Thirty-nine selected candidate TF genes were cloned into the TRV-GFP-PhPDS construct between BamHI and KpnI sites using the Seamless Cloning and Assembly Kit (Novoprotein, Shanghai, China). Oligonucleotide primers used for gene fragment amplification and expression analysis are listed in Supplementary Table S1. The resulting TRV plasmids were electrotransformed into A. tumefaciens strain GV3101. The agro-inoculation was performed according to a previously described protocol (Sun et al., 2017). Three-week-old seedlings of petunia were used for the agro-inoculation. For expression assessment of PhPDS and TF genes, the sequences beyond the inserted fragments were used for primer design (Supplementary Table S1). For GFP fluorescence detection, fluorescent foci were visualized using a portable blue light-emitting diode flashlight as an excitation source at 450 nm (LUYOR-3260RB, LUYOR, Shanghai, China). The images were taken using a digital camera (EOS 20D, Canon, Tokyo, Japan) with a bypass emission filter (LUV-495, LUYOR, Shanghai, China). The intensity of fluorescence was quantified using ImageJ software (NIH, Bethesda, MD, United States).

Agrobacterium-mediated transient assay was performed as previously described (Mußmann et al., 2011) with some minor changes. The complete coding region of PhCOL4 (Supplementary Figure S1) was cloned into pCAMBIA2300 vector between EcoRI and SacI sites. The resulting plasmids were transformed into A. tumefaciens strain LBA4404. The transformed bacteria were cultured in liquid LB media supplemented with 40 mg⋅L^–1 kanamycin for overnight at 28°C. The Agrobacterium cells were harvested and resuspended to an OD600 of 4.0 in the infiltration buffer containing 10 mM MgCl₂, 10 mM MES, and 200 μM acetosyringone. The suspension was infiltrated into fully developed petunia leaves through a syringe. At 6 days post infiltration, the leaves were collected for gene expression analysis.

DAP-Seq binding experiment was carried out according to a previously described method (Yao et al., 2020). In brief, genomic DNA (gDNA) was extracted from petunia leaves using CTAB method, and fragmented through a Covaris M220 (Woburn, MA, United States). Purification of fragmented gDNA was performed using AMPure XP beads (Beckman Coulter, Indianapolis, IN, United States). DNA library was prepared using the NEXTFLEX Rapid DNA-Seq Kit (PerkinElmer, Austin, TX, United States). The coding sequence of PhCOL4 was introduced into pFN19K HaloTag T7 SP6 Flexi expression vector. The TNT SP6 Coupled Wheat Germ Extract System (Promega, Madison, WI, United States) was used to express the Halo-PhCOL4 fusion proteins, which were then captured by Magne Halo Tag Beads (Promega, Madison, WI, United States). The proteins bound to beads were incubated with gDNA fragments at 25°C for 2 h. A high temperature of 98°C for 10 min was used to denature the proteins and release the bound gDNA fragments. The fragments were PCR-amplified, purified, and sequenced. For data processing, all reads were mapped to the reference genome of P. axillaris⁵. Peak analysis was conducted using Macs2 tool. Distribution of peaks in petunia genome sequence was analyzed using Homer (Heinz et al., 2010). Functional annotation of corresponding genes was determined through a blast against NR, NT, Swissprot, and Pfam databases. Protein-bound motifs were predicted using MEME-Chip suite 5.0.5 (Machanick and Bailey, 2011).

Dual luciferase assay was performed using a method as previously described (Sun et al., 2019). The coding sequence of PhCOL4 was ligated into pGreenII 62-SK vector, which was regarded as the effector. To construct the reporters, the promoter sequences of PhRDR6 and PhAGO4 (Supplementary Figure S1) were amplified using the primers listed in Supplementary Table S1, and cloned into pGreenII 0800-LUC vector expressing firefly luciferase (LUC) and Renilla luciferase (REN). The resulting plasmids were transformed into A. tumefaciens strain GV3101. Co-infiltration with effector and reporter constructs was carried out on petunia young leaves. LUC and REN activities were measured using a Tecan Infinite M200 luminometer (Männedorf, Switzerland). The data were expressed as a ratio of LUC to REN.

The 35-bp fragments containing predicted PhCOL4-bound motifs in the PhRDR6 and PhAGO4 promoters were amplified using the primers (Supplementary Table S1). The fragments were ligated into pHIS2 vector to generate the baits. The full-length open reading frame region of PhCOL4 was fused to the GAL4 activation domain in pGADT7-Rec vector as the pray. All constructs were co-transformed into yeast strain Y187, which was grown and spotted on SD-Ura-His-Leu medium supplemented with 100 mM 3-aminotriazole (3-AT) at different dilutions.

To illustrate transcriptional regulatory mechanism of antiviral RNA silencing in petunia, we used an integrative strategy of TRV (PPK20)-infected Illumina RNA-seq and TRV-GFP-PhPDS-based VIGS to identify the TFs regulating the RNA silencing machinery. An important candidate TF was selected for further functional characterization. VIGS and transient overexpression assays were performed to analyze the impact of the identified TF on the expression of RNA silencing components. The target genes downstream of the TF were determined through DAP-Seq, dual luciferase, and yeast one-hybrid assays (Figure 1). Young leaves of P. hybrida at six-leaf-stage were inoculated with TRV (PPK20). A clear symptom development was observed in the inoculated leaves, which exhibited slight leaf yellowing at 3 days post inoculation (dpi) and severe leaf mottling or chlorosis at 6 dpi (Figure 2A). Two regions from TRV (PPK20) genomic RNA1 and RNA2 were chosen for quantification of virus RNA levels (Figure 2B). Virus accumulation analysis showed that virus RNA levels of both TRV-RNA1 and TRV-RNA2 increased gradually through 3 dpi and escalated sharply from 4 to 5 dpi. Then, the levels of TRV RNA transcripts continued to increase but at a much lower speed from 5 to 6 dpi (Figure 2C). Based on the trend of symptom development and viral growth, the infection process of TRV (PPK20) in petunia leaves was divided into three stages: 0 dpi (S0), 3 dpi (S3), and 6 dpi (S6), which were selected for transcriptome analysis.

After data filtration, a total of 141,740,216, 118,866,564, and 117,466,686 clean RNA sequencing reads were obtained for petunia leaves at S0, S3, and S6, respectively. At each stage, 74.2% (105,214,806), 73.1% (86,903,750), and 73.3% (86,084,852) of these reads could be mapped to petunia reference genome. De novo assembly led to the generation of 34,486 unigenes, including 32,928 mapped and 1,558 novel genes (Table 1). The functions of 24,462 mapped unigenes were annotated based on protein databases (Table 1 and Supplementary Table S2). The sequences of all novel transcripts are shown in Supplementary Figure S2.

Compared to S0, S3 displayed 3,413 up-regulated and 1,275 down-regulated unigenes. A large number of DEGs were detected for S6, which has 4,153 up-regulated and 1,782 down-regulated unigenes. When compared to S3, S6 only has a moderate number of DEGs: 895 up-regulated and 784 down-regulated (Figure 2D and Supplementary Tables S3–S5). The DEGs (7,402 in total) from pairwise comparison among S0, S3, and S6 were shown in Figure 2E. The highest number (3,449) of overlapped DEGs was observed between S0-S3 and S0-S6. Consistently, less shared DEGs (745) were found between S0-S3 and S3-S6, and 1,127 overlapped DEGs between S0–S6 and S3–S6. These data suggest that, compared to the control S0, the infected S3 and S6 have relatively closer transcriptome profiles with each other.

Gene Ontology (GO) enrichment analysis was carried out to investigate the functions of identified DEGs. At the molecular function level, the most significantly enriched GO terms correspond to “metabolic process,” “binding,” and “membrane” (Supplementary Figure S3A). At the biochemical pathway level, 1,581 DEGs for S0-S3, 2,007 DEGs for S0-S6, and 428 DEGs for S3-S6 were categorized into 126, 127, and 101 biochemical pathways, respectively. Among these pathways, “protein processing in endoplasmic reticulum,” “Plant–pathogen interaction,” and “Plant hormone signal transduction” were significantly enriched (Supplementary Figure S3B and Supplementary Table S6).

Plant hormones are known to play an important role in defense against viruses. “Plant hormone signal transduction” was detected as significantly enriched pathway in our identified DEGs after TRV (PPK20) infection. Therefore, we measured the levels of common plant hormones ET, ABA, JA, and SA during the TRV (PPK20)-infection process. Results showed that the accumulation of these hormones was significantly boosted by TRV (PPK20) infection. The highest levels of ET and JA in petunia leaves were found at 5 dpi, while ABA and SA contents peaked at 6 dpi (Figure 3A). These results are consistent with our transcriptome data, which showed that ET biosynthetic genes, including PhACO1, PhACO2, PhACO3, PhACO4, PhASC2, PhACS6, and PhACS7, were significantly up-regulated during TRV (PPK20) infection. The transcription of PhZEP, PhNCED3, PhNCED5, PhAAO4, PhSDR2a, PhSDR2b, and PhSDR3b rather than PhNCED7 in the ABA biosynthesis pathway was induced by TRV (PPK20) infection. In addition, most genes involved in JA production exhibited elevated expression except PhLOX3. Several SA biosynthetic genes, such as PhCM1, PhADTs (1 and 6), and PhPALs (1, 2a, and 2b), were also up-regulated upon TRV (PPK20) infection (Figure 3B).

Considering the essential roles of RDR, DCL, and AGO in antiviral RNA silencing, their transcripts differentially expressed during TRV (PPK20) infection were selected from the transcriptome data. Results showed that PhRDR1, PhRDR2, PhRDR6, PhDCL1a, PhDCL1b, PhDCL2, PhDCL3, PhDCL4, PhAGO1a, PhAGO1b, PhAGO2, PhAGO4, PhAGO5, PhAGO8, and PhAGO9 were dramatically up-regulated in TRV (PPK20)-infected petunia leaves. The highest transcriptional changes correspond to PhDCL1b whose expression rose by 130–260 times, followed by a 34–62 times increase in PhDCL2 transcription. In addition, the transcription of PhDCL1a and PhRDR6 increased moderately by less than three times (Figure 4).

To dissect the transcriptional regulatory mechanism for TRV defense in petunia, 469 DEGs encoding putative TFs were identified (Supplementary Tables S7–S9). Compared to S0, there are 286 (204 up-regulated and 82 down-regulated) TFs belonging to 32 families that were differentially expressed at S3. When S6 and S0 were compared, 374 (225 up-regulated and 149 down-regulated) TF-encoding (33 families) DEGs were identified. For S3 and S6, 148 (58 up-regulated and 90 down-regulated) TF-encoding DEGs from 24 families were found. The highest number of TFs included those from bHLH, bZIP, ZFP, ERF, MYB, NAC, and WRKY families (Figure 5A).

To identify critical transcriptional regulators involved in antiviral RNA silencing, 39 TRV-induced TFs were selected for TRV-GFP-PhPDS-based VIGS assay (Table 2). The fragments of these TFs were cloned into the multiple cloning sites of TRV-GFP-PhPDS vector to generate the recombinant constructs (Figure 5B). GFP and PhPDS were used as reporter genes for virus accumulation and gene silencing. Normal green fluorescence and leaf photobleaching were observed due to GFP expression and PhPDS silencing. Compared to TRV-GFP-PhPDS control, the inoculation with TRV-GFP-PhPDS constructs containing the fragments of PhARF11, PhbHLH41, PhbHLH93, PhZPT4-3, PhCOL4, PhGATA11, PhHSF-B3A, PhLOB1, PhNAC90, and PhWRKY75 resulted in enhanced fluorescent signals at 4 dpi. In contrast, the fluorescent brightness was suppressed in the leaves infected with TRV-GFP-PhPDS constructs containing PhPIF3, PhZAT10, PhERF22, PhERF1B, PhHSF-B3B, PhHSF-A2, and PhWRKY40 inserts (Figure 6). In addition, TRV-GFP-PhPDS/TGA2.1-, ZFP6-, and ERF22-infected systemic leaves showed promoted leaf photobleaching at 14 dpi, in comparison with the control. However, the plants inoculated with TRV-GFP-PhPDS vectors containing PhbHLH41, PhbHLH93, PhZPT4-3, PhZAT10, PhCOL4, PhBBX20, PhERF3, PhKD1-3, PhHSF-B3A, PhHSF-B3B, PhNAC90, PhJUB1B, and PhWRKY75 fragments displayed impaired photobleaching phenotypes (Figure 6).

Transcript levels of 39 candidate TFs were examined in the uppermost leaves of petunia plants infected with TRV constructs. Results showed that all TF transcripts were significantly down-regulated in the corresponding TRV construct-infected petunia leaves (Supplementary Figure S4). Compared to TRV-GFP-PhPDS control, higher levels of GFP fluorescence intensity were observed with the down-regulation of PhARF11, PhbHLH41, PhbHLH93, PhZPT4-3, PhCOL4, PhGATA11, PhHSF-B3A, PhLOB1, PhNAC90, and PhWRKY75. On the contrary, silencing of PhPIF3, PhZAT10, PhERF22, PhHSF-B3B, PhHSF-A2, and PhWRKY40 led to lower fluorescence levels. TRV-RNA1 and TRV-RNA2 accumulation levels exhibited highly similar patterns to GFP fluorescence intensities (Figure 7A).

At the phenotype level, TRV-GFP-PhPDS/TGA2.1-, ZFP6-, and ERF22-infected plants displayed remarkably higher rates of photobleached leaves than the control (TRV-GFP-PhPDS). In contrast, TRV-GFP-PhPDS/bHLH41-, bHLH93-, ZPT4-3-, ZAT10-, COL4-, BBX20-, ERF3-, KD1-3-, HSF-B3A-, HSF-B3B-, NAC90-, JUB1B-, and WRKY75-infiltrated plants demonstrated much lower rates of photobleached leaves. Compared to the control, the transcription of PhPDS was lower in the leaves with enhanced photobleaching, but higher in the leaves with reduced photobleaching (Figure 7B), which is consistent with the leaf phenotype.

In view of the prominent impact of PhCOL4 (SGN accession no. Peaxi162Scf00045g01824) silencing on GFP fluorescence and leaf photobleaching in the VIGS assay, we selected it for further functional characterization. Expression levels of PhCOL4 were examined in its silenced and transiently overexpressed petunia leaves. At 4 dpi with TRV-GFP-PhPDS/COL4, petunia leaves exhibited decreased transcript levels of PhCOL4 and several RNA silencing-related genes, including PhRDR6, PhDCL2, PhDCL4, PhAGO2, and PhAGO4, compared to TRV-GFP-PhPDS control (Figure 8A). Transient overexpression of PhCOL4 in petunia leaves resulted in increased transcript abundances of PhRDR6, PhDCL2, PhDCL4, PhAGO2, and PhAGO4, compared to empty vector control (Figure 8B).

To screen downstream target genes of PhCOL4, a DAP-Seq assay was performed. A total of 94,863,846 high-quality reads were generated from PhCOL4-bound products, compared to input control (Figure 9A). A number of peaks were located before transcription start sites of genes (Figure 9B). Of the 1,626 peaks, 660 (40.59%) were located within the 2kb upstream promoter regions. The remaining peaks were located within the regions of 5′UTR, 3′UTR, exon, intron, and downstream (2 kb after transcription end sites) (Figure 9C). Analysis of MEME and DREME revealed two top PhCOL4-bound motifs with the core element TTCTT (or AAGAA in the opposite strand) (Figure 9D). Through the functional annotation of peaks, two important genes, PhRDR6 and PhAGO4, involved in the RNA silencing pathway were found (Figure 9E).

To determine the transactivation of PhRDR6 and PhAGO4 promoters by PhCOL4, a dual luciferase assay based on the effector and reporter constructs was carried out (Figure 10A). Compared to empty vector control, the co-expression of 35S:PhCOL4 with pPhRDR6:LUC or pPhAGO4:LUC caused a dramatic increase in LUC activity (Figure 10B). To further confirm the interaction of PhCOL4 with predicted motifs by DAP-Seq, a 35-bp fragment containing motif sequence in the PhRDR6 or PhAGO4 promoter was used as a probe for yeast one-hybrid assay. The constructs of promoter and PhCOL4 were used as the bait and pray, respectively (Figure 10C). The yeast cells co-transformed with PhCOL4 and wild-type promoter fragments grew well on the selective medium in the presence of 100 mM 3-AT. However, the growth of yeast transformants with mutant promoter fragments were remarkably inhibited (Figure 10D). These results revealed the direct interaction between PhCOL4 and the promoters of PhRDR6 and PhAGO4.