The vectors pDONR207, pEAQ-HT-DEST3, and Agrobacterium tumefaciens C58C1 were provided by Professor Juan Antonio García from the National Centre for Biotechnology (CNB), Spain. Escherichia coli DH5α competent cells were purchased from TaKaRa Bio (Osaka, Japan). The infection clone of TVMV was described by Zhao et al. (2020).

Seeds of N. benthamiana plants were sown in soil, and the plants were cultured in a greenhouse at 22°C under a 16/8 h light/dark cycle. In general, 4 ~ 5-leaf-stage N. benthamiana plants with uniform growth and size were selected for the experiments.

Specific primers were designed according to the sequences of the 3'UTR of TVMV, GFP, and the expression vector pEAQ-HT-D3. The primers used are listed in Supplementary Table 1. High-fidelity PrimeSTAR GXL DNA polymerase (TaKaRa) was used to amplify fragments in overlap PCR. The TVMV infection clone, pLXB-TVMV, was described by Zhao et al. (2020). To facilitate fragment recombination, an overlapping segment of 24 bp in length was designed at the 3' and 5' ends of each primer. Two primer pairs (attB-3'UTR-F1/3UTR-GFP-R1 and 3'UTR-GFP-F2/attB-GFP-R2) were used to amplify the PCR products of the TVMV 3'UTR and GFP using pLXB-TVMV and pEAQ-GFP (Sun et al., 2023) as templates. The PCR fragments of the TVMV 3'UTR and GFP were mixed at a ratio of 1:1 as the template, and primers attB-3'UTR-F1 and attB-GFP-R2 were used to amplify the3'UTR-GFP fragment using TaKaRa LA Taq DNA polymerase. The three PCR products were detected and purified using 1% agarose gel electrophoresis.

The expression vector pEAQ-3'UTR-GFP was constructed by Gateway technology. The obtained PCR product of 3'UTR-GFP was subjected to BP reaction to construct an entry clone. Then, 3.5 μl purified PCR product, 0.5 μl donor vector pDONR207, and 1 μl BP enzyme mix were fully mixed and incubated at 25°C for 6 h. Subsequently, 5 μl of all BP reaction products was transformed into E. coli DH5α. The positive antibody was screened on LB solid medium containing 50 μg/ml gentamicin. Plasmid pDONR-3'UTR-GFP was extracted and detected with restriction endonucleases Eag I and Pst I. The plasmid was purified and sequenced using BGI Group (Shenzhen, China).

The obtained pDONR-3'UTR-GFP and the target vector PEAQ-HT-DEST3 were used for LR reaction. Then, 0.35 μl pDONR-3'UTR-GFP (100 ng), 1 μl target vector pEAQ-HT-DEST3 (100 ng), 1 μl LR enzyme mix, and 2.65 μl ddH₂O were fully mixed and incubated at 25°C for 6 h. Next, 5 μl of all LR reaction products was transformed into E. coli DH5α. The positive antibody was screened on LB solid medium containing 50 μg/ml kanamycin. The plasmid PEAQ-3'UTR-GFP was extracted, purified, and sequenced by BGI Group.

The constructed expression vector, pEAQ-3'UTR-GFP, was transformed into Agrobacterium C58C1 and expressed transiently in N. benthamiana. GFP expression was detected in the leaves. The cultured Agrobacterium C58C1 cells containing pEAQ-3'UTR-GFP were suspended in induction buffer (0.01 M MES solution, 0.01 M MgCl₂, 0.15 mM Aceto mixed solution). The induced pEAQ-3'UTR-GFP bacterial suspension was injected into the third and fourth leaves of N. benthamiana from the back of the leaves with a 1 ml single-use aseptic injector to fully infiltrate the leaves. Then, the inoculated N. benthamiana was cultured at 22°C. N. benthamiana injected with pEAQ-GFP and pEAQ-HT-DEST3 were used as controls, with six plants per treatment. The expression of GFP in the leaves was observed 48 h after infection. The inoculated leaves were collected and flash frozen in liquid nitrogen to detect GFP expression.

Total protein was extracted using liquid nitrogen urea lysate via high-speed centrifugation. The plant samples were ground in a precooled mortar and pestle and collected in 1.5 ml centrifuge tubes. Protein lysate (pH = 7.5; 125 mM TrisHCl, 20% sodium dodecyl sulfate (SDS), 6 M urea, 5% β-mercaptoethanol, and a small amount of bromophenol blue) was added to the centrifuge tube at a 2:1 ratio, and the mixture was incubated in a 95°C metal bath for 10 min, followed by 2 min on ice. The centrifuge tubes were then centrifuged at 13,400 × g for 10 min at 4°C. Five microliters of the prepared protein samples were subjected to sodium dodecyl sulfate–polyacrylamide gel electrophoresis. Following electrophoresis, the membrane transfer process was initiated, and the protein was transferred to a nitrocellulose membrane and blocked in a blocking solution for 2 h. The nitrocellulose membrane was washed with TBST and incubated with 5,000-fold-diluted GFP monoclonal antibody (Bio-Rad, Hercules, California, USA), before being washed again with TBST and incubated with 5,000-fold-diluted alkaline phosphatase-labeled goat anti-rabbit IgG (Bio-Rad, Hercules, California, USA) secondary antibody. Then, the membrane was washed a third time with TBST and treated with an ECL color development solution (Bio-Rad, Hercules, California, USA) to visualize the protein bands on an Odyssey® Fc near-infrared dual-color laser and chemiluminescence dual-imaging function imaging system.

TRIzol reagent was used to extract total RNA and further synthesize cDNA (TransGen Biotech, Beijing, China). Based on the sequence of the GFP coding region, Primar5.0 was used to design qPCR primers. The upstream primer was GFP (el)-F, and the downstream primer was GFP (el)-R. GFP expression was detected by quantitative real-time PCR (RT-qPCR). NtUB1 was used as the internal reference gene, and its upstream and downstream primers were NtUB1-F and NtUB1-R, respectively (Supplementary Table 1). Three replicates were analyzed for each gene. The 2^−ΔΔCT method was used to calculate the level of gene expression.

Total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific, Waltham, Massachusetts, USA) following the manufacturer’s instructions. Total RNA quantity and purity were determined using a Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent Technologies, Santa Clara, California, USA), and high-quality RNA samples with RNA integrity number > 7.0 were used to construct the sequencing library. RNA libraries were sequenced on an Illumina Novaseq 6000 platform in PE150 mode by LC Bio Technology Co., Ltd. (Hangzhou, China). The transcriptome comprised nine cDNA libraries, consisting of three biological replicates and three sets of treatments.

Raw data generated by sequencing were preprocessed and filtered to obtain clean data. Clean data were compared with the N. benthamiana genome¹ to obtain comprehensive alignment information. Simultaneously, according to the gene location information specified in the genome annotation file gtf, the following statistics were performed: read comparison between the sequencing data and the reference genome and the regional distribution analysis of the sequencing data and the reference genome. The quantitative status of each sample gene was determined using StringTie. After quantitative gene analysis, differentially expressed genes (DEGs) were screened for functional annotation and enrichment analysis. Gene expression was expressed as fragments per kilobase million. Fold change ≥2 and false discovery rate < 0.01 were used as screening criteria. The raw sequence data were submitted to the NCBI Short Read Archive (SRA) with accession number PRJNA1134345.

EdgeR was used to analyze the DEGs in N. benthamiana after the transient expression of 3'UTR-GFP, GFP, and Vector for 48 h. The number of up-and downregulated genes was obtained by using |log₂FC| ≥ 1 and q < 0.05 as the screening criteria. Comparison and annotation were performed using the Gene Ontology (GO) database, and GO functional annotations were obtained using ‘Blast2GO’ software (Conesa et al., 2005). All GO functional annotations were classified using ‘WEGO’ software (Ye, 2006), and the DEGs were assigned to cellular component, biological process, and molecular function categories in the GO database. Pathway enrichment analysis was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, the main public pathway database, which primarily analyzes the metabolic pathways and functions of cellular proteins. The KEGG database was used to annotate the metabolic pathways of each independent gene translational protein.

Nine DEGs of N. benthamiana after transient 3'UTR-GFP expression were selected for RT-qPCR verification, and the necessary primers were synthesized by Shanghai Sangon Biological Co. (Supplementary Table 1). The NtUB1 gene was used as an internal reference gene, and cDNA templates were obtained from the transcriptome sequencing samples. The 2^−ΔΔCT method was used to determine the relative expression.

The 3'UTR-GFP construction was conducted in accordance with the Gateway system (Figure 1A). To obtain the 3'UTR-GFP fragment, fragments of 3'UTR (352 bp) and GFP (772 bp) of TVMV were amplified (Figures 1B,C). Then, the 3'UTR of TVMV was fused upstream of GFP by overlap PCR to create the 3'UTR-GFP (1,034 bp; Figure 1D). Next, 3'UTR-GFP was introduced into the expression vector pEAQ-HT-D3 via the Gateway recombination system. First, the 3'UTR-GFP fragment was cloned into the entry clone pDONR207 by BP reaction. pDONR-3'UTR-GFP was digested by Eag I and Pst I, which resulted in three bands of 4,301, 2,575, and 1,726 bp in length (Figure 1E). Then, the 3'UTR-GFP fragment was introduced into the expression vector pEAQ-HT-D3 via LR reaction, resulting in the creation of pEAQ-3'UTR-GFP. The plasmids pDONR-3'UTR-GFP and pEAQ-3'UTR-GFP were subjected to sequencing, and the 3'UTR-GFP sequences in each plasmid were subsequently analyzed (Supplementary Figure 1; Supplementary Table 2). The 3'UTR-GFP sequence was found to be accurate and uniform across all plasmids, as evidenced by the sequencing peak figures (Supplementary Figure 2; Supplementary Table 2). In order to check if the 3'UTR structure in pEAQ vector is the same as that of TVMV virus. We have predicted the structure of 3'UTR in vitro by RNAfold Webserver and make a comparison with the structure in TVMV viral RNA. We have used the 253 nt of 3'UTR to perform the structure prediction. Regarding with the TVMV virus, the limited length of RNA for structure prediction is 7,500 nt. We choose the 7,500 nt of TVMV calculated from 3' end. It appears that the structure of 3'UTR alone and that of in virus are similar (Figure 2). It is supposed that the 253 nt of 3'UTR will be transcribed in pEAQ-D3-3'UTR after agro-infiltrated in N. benthamiana.

The pEAQ-3'UTR-GFP expression vector was transformed into Agrobacterium C58C1 and inoculated into leaves of N. benthamiana. GFP expression was observed 48 h after injection. Pronounced green fluorescence was evident on N. benthamiana leaves that were injected with pEAQ-GFP, whereas the fluorescence on leaves injected with pEAQ-3'UTR-GFP was considerably weaker (Figure 3A).

The inoculated leaves of N. benthamiana were harvested, and the protein extracts were utilized directly for Western blotting with anti-GFP antibody. The results demonstrated the absence of GFP protein in the injected leaf samples inoculated with the empty vector pEAQ-HT-D3. The level of GFP expression in samples treated with pEAQ-3'UTR-GFP was lower than that in samples receiving pEAQ-GFP (Figure 3B). The GFP signal was quantitatively analyzed, and the results are shown in Figure 3C. The GFP signal in the pEAQ-3'UTR-GFP samples exhibited a markedly diminished intensity in comparison to that observed in the pEAQ-GFP samples (p = 0.004; Figure 3C).

RT-qPCR was performed to confirm the difference in GFP expression between the two samples. GFP expression was in accordance with the results of Western blotting, as shown in Figure 3D. The GFP expression in pEAQ-3'UTR-GFP samples was significantly lower than that in pEAQ-GFP samples (p = 0.006).

Transcriptome sequencing analysis of N. benthamiana leaves infiltrated with empty vector pEAQ-HT-D3, pEAQ-3'UTR-GFP, and pEAQ-GFP was performed using an Illumina Novaseq™ 6000 platform (LC Bio Technology CO., Ltd. Hangzhou, China). The three treatments yielded 5.17, 5.36, and 6.19 G clean reads, respectively. The G + C content of the bases was >42.5%, and the proportion of Q30 bases was >98.20% (Supplementary Table 3). The genome sequence of N. benthamiana was used as a reference. The clean reads of all the samples showed 89.63–92.39% similarity with the reference genome, as shown in Supplementary Table 3. These findings confirm the accuracy and reliability of the sequencing assembly results, making them suitable for further analyses.

The distribution of gene expression values across the samples revealed that the nine boxes were similar in size, indicating a high degree of similarity in gene density among the nine samples. The consistency in whisker lengths demonstrated that samples exhibited similar variance. The median line of the samples was closer to the high values of the boxes, with a long whisker shape extending towards the low values, thereby producing a negative bias (Supplementary Figure 3A). A correlation analysis of the gene expression information from the 3'UTR-GFP, GFP, and Vector samples revealed that the Pearson correlation coefficient for all three groups of samples was above 0.9 (Supplementary Figure 3B), indicating a high degree of correlation in the expression patterns between the three groups. This confirms the reliability of the data, which can be used for further analyses.

DESeq2 (Anders and Huber, 2010) was used to analyze the DEGs among 3'UTR-GFP vs. GFP, 3'UTR-GFP vs. Vector, and GFP vs. Vector samples, revealing 159, 777, and 352 DEGs, respectively (|log₂FC| ≥ 1 and q < 0.05; Figure 4A; Supplementary Table 4). Overall, 66, 356, and 232 DEGs were upregulated and 93, 421, and 120 were downregulated in the 3'UTR-GFP vs. GFP, 3'UTR-GFP vs. Vector, and GFP vs. Vector comparisons, respectively (Figures 4C–E). Cluster analysis identified 10 DEGs common to all comparisons and 50, 478, and 120 unique DEGs in the 3'UTR-GFP vs. GFP, 3'UTR-GFP vs. Vector, and GFP vs. Vector comparisons (Figure 4B; Supplementary Table 5).

The 159 DEGs in the 3'UTR-GFP vs. GFP comparison were subjected to GO functional enrichment analysis and quantitative statistics (Figure 5A). We found 25 categories of biological processes, all of which showed differential gene enrichment. Among them, “regulation of transcription, DNA-templated” genes were significantly enriched, with one gene was upregulated in 3'UTR-GFP group relative to its expression in GFP group. Fourteen illogical process genes were downregulated in 3'UTR-GFP group relative to its expression in GFP group, accounting for 9.4% of the total downregulated genes. We identified 22 categories of cellular components, 10 of which showed differential gene enrichment. The “nucleus” and the “integral component of membrane “had the most gene enrichment, with one upregulated gene each and two and one downregulated genes, respectively. Differential gene enrichment was also observed for all 22 molecular function classifications. In particular, genes in “DNA binding,” “DNA-binding transcription factor activity,” “ATP binding,” and “catalytic activity” had the most gene enrichment, with one, one, six, and three upregulated genes, respectively. These genes accounted for 1.3, 1.3, 8, and 4% of the total upregulated genes, respectively. The number of downregulated genes was 16, 13, 6, and 8, accounting for 10.7, 8.7, 4.0, and 5.4% of the total downregulated genes, respectively.

The 777 DEGs in the 3'UTR-GFP vs. Vector comparison were also subjected to GO function enrichment analysis and quantitative statistics (Figure 5B). We found 25 biological process categories showing differential gene enrichment. Among them, “regulation of transcription, DNA-templated” genes had the most gene enrichment, with 27 upregulated and 45 downregulated genes in the 3'UTR-GFP group, accounting for 4.4 and 6.3% of the total up-and downregulated genes, respectively. All 22 cell components categories showed differential gene enrichment. Membrane gene enrichment was significant; the number of upregulated and downregulated genes was 15 and 17, respectively. Differential gene enrichment was also observed for all 22 molecular function classifications, with significant enrichment in genes related to “DNA binding,” “protein binding,” “DNA-binding transcription factor activity,” and “ATP binding,” with 23, 40, 20, and 31 upregulated genes accounting for 3.8, 6.5, 3.3, and 5.1% of the total upregulated genes, respectively. The number of downregulated genes was 46, 23, 41, and 13, representing 6.4, 3.2, 5.7, and 1.8% of the total downregulated genes, respectively.

Finally, GO functional enrichment analysis and quantitative statistics were performed for the 352 DEGs in the GFP vs. Vector comparison (Figure 5C). Transcriptional regulation and DNA template genes had the most gene enrichment, with 23 upregulated genes accounting for 5.0% of the total upregulated genes. Eleven genes were downregulated, representing 6.6% of all the downregulated genes. Among the 22 cell component categories, 16 were enriched in DEGs. Membrane genes had the most gene enrichment, with 10 upregulated and 15 downregulated genes. Differential gene enrichment was observed for all 22 molecular function classifications. Genes related to “DNA binding,” “ATP binding,” “DNA-binding transcription factor activity,” and “protein binding” had the most gene enrichment, with 10, 8, 6, and 4 upregulated genes accounting for 2.2, 1.8, 1.3, and 0.9% of the total upregulated genes, respectively. In contrast, five, two, zero, and one genes were downregulated, respectively.

KEGG pathway analysis of the DEGs in the 3'UTR-GFP group vs. the GFP group annotated 98 DEGs (61% of all genes), of which 91 had the most gene enrichment. The enrichment level of DEGs in the plant–pathogen interaction pathway was the highest, with 18 DEGs enriched, 3 of which were upregulated and 15 downregulated (Figure 5D; Supplementary Table 6).

By removing the DEGs in the GFP vs. Vector comparison from those obtained in the 3'UTR-GFP vs. GFP and 3'UTR-GFP vs. Vector comparisons, we identified 611 DEGs related to 3'UTR response (Supplementary Table 7). These 611 DEGs were subjected to GO functional enrichment analysis, with results presented in Figure 6A. Overall, 350 genes were significantly enriched in 41 GO terms (p < 0.05). In the biological process branch, the two terms containing the most DEGs were “transcriptional regulation, DNA-templated” and “carbohydrate metabolic processes.” In the cellular component branch, the three terms containing DEGs were “cell wall,” “apoplast,” and “anchored component of the membrane.” In the molecular function branch, the terms containing the most DEGs were “DNA binding,” “sequence-specific DNA-binding transcription factor activity,” and “heme-binding.”

The results of KEGG signaling pathway enrichment showed that 269 DEGs were enriched in 95 KEGG pathways, of which 6 were significantly enriched (p ≤ 0.05). As shown in Figure 6B, the plant–pathogen interaction (ko04626) pathway had the highest enrichment level of DEGs, with 44 DEGs annotated to the pathway, of which 11 were upregulated and 33 were downregulated in the 3'UTR-GFP group (Figure 6C; Supplementary Table 8).

Among the 611 DEGs related to the 3'UTR reaction, we obtained 80 related to plant resistance, of which 29 were upregulated and 51 were downregulated (Figure 7A; Supplementary Table 9). The upregulated genes included those encoding defensin-like protein 3, acidic endochitinase, peroxidase, respiratory burst oxidase homolog protein E, disease resistance protein, WRKY transcription factor 6, calcium-binding protein, calmodulin-binding protein, ethylene-responsive transcription factor 4, protein kinase, protein phosphatase 2C 35, potassium channel subfamily K member 3, and ferric reduction oxidase 5. The downregulated genes included those encoding endoglucanase, beta-glucan-binding protein 4, peroxidase, WRKY transcription factor, calcium-binding family protein, ethylene-responsive transcription factor, protein kinase, protein phosphatase 2C family protein, kinesin-related protein 1, Myb-related protein 306, and histone H2A.1.

Nine DEGs, namely, peroxidase 5 (POX), WRKY transcription factor 6 (WRKY6), receptor kinase 2 (RK), acidic endochitinase (AE), calcium-binding EF-hand family protein (CbFP), ethylene-responsive transcription factor 10 (ERF), histone H2A.1 (H2A), protein phosphatase 2C family protein (PP2C), and WRKY transcription factor 22 (WRKY22), were selected for RT-qPCR verification. Each gene signal was quantitatively analyzed, and the results are shown in Figure 7B. Compared with their expression in the Vector group, four genes were upregulated and five were downregulated in the 3'UTR-GFP group, of which three upregulated genes were significantly differentially regulated. The Log₂FC value was calculated using the 2^−ΔΔCT method. The results showed that the expression patterns of the nine selected genes were consistent with the transcriptome data compared to the RNA-seq data (Figure 7C). This finding corroborated the reliability of the transcriptome data.