TSWV (isolate: GC-1) was isolated from lettuce (Lactuca sativa) in Yunnan Chenggong, and HCRV (isolate: HLS1-2) was isolated from susceptible spider lily (Hymenocallis littoralis) in Kunming. TSWV or HCRV was inoculated and maintained in Emilia sonchifolia (L.) plants. Six- to eight-week-old plants were used for viral inoculation. TSWV or HCRV-infected E. sonchifolia leaves were collected and ground in an inoculation buffer to inoculate N. benthamiana plants. Tobacco (N. benthamiana) seedlings with four to five true leaves were divided into four groups, which included inoculation with TSWV, HCRV, and TSWV + HCRV. Buffer inoculation was used as the control group (CK). The inoculum was prepared by grinding 1 g of systemically infected fresh leaves in 10 ml of inoculation buffer (Ocampo et al., 2016). The identified symptomatic young leaves were collected and ground in a sterilized mortar and then dipped in cotton swabs on the tobacco leaves for inoculation. The symptoms were recorded after inoculation at 7 days post-inoculation (dpi). The leaves above the inoculated leaves were collected at 1, 7, and 14 dpi and stored at −80°C.

Total RNA from different treatments of tobacco with three biological replicates was extracted using a TRIzol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The rRNAs were then removed to retain the mRNAs and non-coding RNAs (ncRNAs). The enriched mRNAs and ncRNAs were fragmented into short fragments using fragmentation buffer and reverse transcribed into cDNA with random primers. Next, the cDNA fragments were purified with a QIAquick PCR Extraction Kit (Qiagen, Venlo, Netherlands), end-repaired, poly(A) added, and ligated to Illumina (San Diego, CA, USA) sequencing adapters. The digested products were sequenced using Illumina HiSeq™ 4000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Total cDNA was prepared from RNA using the Superscript III™ Reverse Transcriptase kit (Invitrogen, Cat. 18080–044). A one-tenth dilution of the reverse transcription final reaction was prepared; 1 μl of the dilution was used as a template for the qPCR consisting of 0.4 μM of each primer and 1× SYBR Green Master Mix (QuantiTect^® SYBR^® Green PCR Kit; Qiagen, Cat. 2041453). Transcript levels were normalized to the expression level of Actin.

The short read alignment tool Bowtie2v2.4.4 (Langmead and Salzberg, 2012) was used to map reads to a ribosomal RNA (rRNA) database. The reads were further filtered by FASTPv0.23.1 (Chen et al., 2018), so that high-quality clean reads were obtained (Chen et al., 2018). An index of the reference genome was then constructed and paired-end clean reads were mapped to the reference genome using HISAT2v2.1.0 (Kim et al., 2015) with “-rna-strandness RF.” The reconstruction of transcripts was performed with the software StringTiev2.2.0 (Pertea et al., 2016) combined with HISAT2. Novel transcripts were then aligned to the nonredundant (Nr), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Ontology (GO) databases to annotate functional proteins.

Transcript abundances were quantified using the software StringTiev2.2.0 in a reference-based approach. For each transcription region, a fragments per kilobase of transcript per million mapped reads (FPKM) value was calculated to quantify the abundance of expression. The differentially expressed transcripts of coding RNAs and long non-coding RNAs (lncRNAs) were analyzed, respectively. Differential expression analysis of the RNAs and lncRNAs was performed using DESeq2v3.11 (Li and Dewey, 2011) software between the two different groups. The genes/transcripts with the parameter of false discovery rate (FDR) below 0.05 and absolute fold change ≥2.0 were considered differentially expressed genes/transcripts.

The N. benthamiana plantlets under normal conditions were used for RNA extraction. The 5’ RLM-RACE adapter was obtained through RNA processing following the instructions in the book. Then, the cDNA acquired by reverse transcription was used as the template for nested PCR according to the protocol of the FirstChoice^® RLM-RACE kit (no. AM1700; Invitrogen). The cloned products were purified for further sequencing.

To construct TRV vectors, 350 bp DNA fragments of target genes of N. benthamiana were cloned using a primer pair containing EcoRI and BamHI sites and then ligated to pTRV2. The correct recombinant vectors were transformed into Agrobacterium GV3101, then injected into tobacco. The empty vector-containing Agrobacterium solution was used as a control. At 14 dpi, the transcript level of target genes in plants was assessed using RT-qPCR. The silenced tobaccos were inoculated with HCRV and TSWV.

A principal component analysis (PCA) was performed with the R package gmodels v2.18.1.1 (http://www.r-project.org/). The correlation coefficient between the two replicas was then calculated to evaluate the repeatability between samples based on the Pearson value.

Gene Ontology and KEGG pathway enrichment analyses were performed using the R package Clusterprofilerv4.0 (Wu et al., 2021). GO has three ontologies, which include molecular function, cellular component, and biological process. Significantly enriched GO or KEGG pathway terms in the DEGs that were compared with the genome background were defined by a hypergeometric test. The p-values were calculated using FDR correction and taking FDR ≤0.05 as a threshold. The GO enrichment results were visualized using the online tool omicshare (https://www.omicshare.com/), while the scatter plot of KEGG pathway results was generated using ggplot2v3.3.0 (Ito and Murphy, 2013).

MapMan software v3.1 (Thimm et al., 2004) was used to visualize cacao gene expression data in the context of metabolic pathways. MapMan uses a plant-specific ontology that classifies genes into well-defined hierarchical categories, denominated BINs. Cacao genes were assigned to BINs using the Mercator automated annotation pipeline (http://mapman.gabipd.org/web/guest/mercator). Differentially represented MapMan pathways were defined by a two-tailed Wilcoxon rank sum test corrected by the Benjamin–Hochberg method (FDR = 0.05).

The protein–protein interaction network was identified using String v10 (Subramanian et al., 2005) and visualized using Cytoscape software v3.7.2 (Shannon, 2003). In the PPI network, nodes represent the target proteins, while edges represent the predicted or validated interactions between proteins.

The topological network of miRNA and its targets was visualized using Cytoscape v3.7.2 (Shannon, 2003). Topological analysis of target genes was performed using the NetworkAnalyzer plug-in and CytoNCA plug-in of Cytoscape. Target proteins were filtered separately according to betweenness centrality (BC), closeness centrality (CC), and degree centrality (DC), which were calculated using the CytoNCA plug-in. The top 10 miRNAs of each subnetwork were retrieved, and overlapped genes were selected as key targets in the present research.

The total RNA molecules in the size range of 18–30 nt were enriched by polyacrylamide gel electrophoresis (PAGE). The 3’ adapters were then added, and the 36–48 nt RNAs were enriched. The 5’ adapters were then ligated to the RNAs as well. The ligation products were reverse transcribed by PCR amplification, and the 140–160 bp size PCR products were enriched to generate a cDNA library and sequenced by Gene Denovo Biotechnology Co. using Illumina HiSeq X. All the clean tags were aligned with small RNAs in the GenBank database (Release 209.0) to identify and remove rRNA, small coding RNA (scRNA), small nucleolar RNA (snoRNA), small nuclear RNA (snRNA), and tRNA. All the clean tags were also aligned with the reference genome.

All the clean tags were then searched against the miRBase database (Release 22) to identify known (species studied) miRNAs. All the unannotated tags were aligned with the reference genome. The level of miRNA expression was calculated and normalized to transcripts per million (TPM). The differential expression of miRNAs was analyzed using edgeR software v4.2 (Robinson et al., 2010) between the two different groups. We identified miRNAs with a fold change of ≥2.0 and P <0.05 in a comparison as significant differentially expressed miRNAs. Based on the sequences of the existing miRNAs, known miRNAs, and novel miRNAs, the candidate target genes were predicted using the software PatMatch (Version 1.2).

Full-length miRNAs (Accession numbers: novel-m0693-5p, novel-m0782-3p, novel-m0915-3p, novel-m1488-5p, and novel-m1992-3p) of N. benthamiana were cloned using specific primers to construct overexpression vectors. The DNA products were double-digested with EcoRI and BamHI and ligated to the pBIN121-MYC vector restricted by the same enzymes. The correct recombinant vectors were confirmed by sequencing.

Positive colonies of Agrobacterium GV3101 harboring overexpression vectors or empty vectors were cultured in 5 ml of LB liquid medium plus 50 mg/L kanamycin and 20 mg/L rifampicin with shaking (200 rpm) at 28°C to an OD 600 value of 2.0, followed by centrifugation at 8,000 rpm for 10 min. The Agrobacterium cells were re-suspended to an OD 600 of 0.6 in MES buffer (10 mM MES-KOH, pH 5.2, 10 mM MgCl₂, and 100 µM acetosyringone). The overexpression of miRNAs in tobacco was performed through agro-infiltration. Then, the expression of resistance genes in tobacco with miRNA overexpression was determined using RT-qPCR, and the experiments were performed three times.

In this study, we found that both TSWV and HCRV could infect the leaves of N. benthamiana, leading to observable lesions on the related leaves at 7 dpi ( Figure 1A ). Symptoms on tobacco leaves began to appear at 15 days post-inoculation (dpi) of individual HCRV, including chlorosis, deformity, and necrosis of leaves ( Figure 1A ). Meanwhile, a single TSWV inoculation could cause symptoms on tobacco leaves at 7 dpi, and it mainly induced concentric ring lines, chlorosis, yellowing, and necrosis of leaves ( Figure 1A ). Notably, TSWV and HCRV could synergistically infect the tobacco leaves as a complex and induce more observable lesions on the leaves at 7 dpi compared with the individual infection of TSWV or HCRV; the biomass of TSWV and HCRV was the highest in tobacco at 7 dpi of co-infection ( Figure 1A ). The incidence rate of lesions caused by co-infection of the TH complex was faster than that of single HCRV or TSWV ( Figure 1A ). TH complex inoculation mainly caused symptoms on tobacco leaves at 5–7 dpi, and it mainly induced symptoms represented by concentric rings, chlorosis, yellowing, and necrosis of leaves ( Figure 1A ).

Subsequently, the transcriptomes of N. benthamiana inoculated with TSWV (T), HCRV (H), and the TSWV–HCRV (TH) complex during different disease response stages (1, 7, and 14 dpi) were sequenced using an Illumina platform. Healthy leaves without infection were used as a control (CK). A total of 264.8 Gb of clean reads were obtained from the 36 samples, which were included in the Illumina sequencing data from all the treatments. After filtering the low-quality reads, polyA, and adaptor, we obtained high-quality clean reads in the transcriptome profiles ( Figure 1B ). All the clean reads from each sample were then mapped uniquely to the N. benthamiana reference database. We compared all the samples to identify the significant differentially expressed genes (DEGs) using the threshold FDR ≤0.05 and log 2 |fold change| >1 ( Figure 1C ). The results showed that all the TSWV, HCRV, and TH complexes induced dramatic variations in the transcription pattern of N. benthamiana ( Figure 1C ). There were 1,574 upregulated and 2,844 downregulated genes in N. benthamiana following 1 dpi of HCRV. TSWV induced the differential expression of more genes in N. benthamiana at 14 dpi ( Figure 1C ). Consistently, the most significant variations in gene expression caused by TSWV infection were also identified at 14 dpi ( Figure 1C ). Additionally, the variations induced by TH complex infection were smaller than those induced by TSWV and HCRV ( Figure 1C ). Although the expression differences between T1, H1, and TH1 were small, the TH samples exhibited significant changes in gene expression at 7 dpi compared with T and H ( Figure 1C ). Thus, we propose that the HT complex could alter the pattern of expression of N. benthamiana to escape the plant defense responses and contribute to their ability to infect.

To investigate the influences of TSWV, HCRV, and TH complex infection on the dynamic transcriptome pattern of N. benthamiana, a principal component analysis (PCA) was performed on all the transcriptome profiles from each treatment with three biological replicates ( Figure 2 ). We identified that the dynamic expression patterns of tobacco following HCRV infection were significantly affected by TSWV infection during its infection stage, and the samples after 1, 7, or 14 dpi of TSWV were distributed in different regions ( Figure 2A ). Principal component 1 explained 43.6% of the differences among all the samples, and these differences were primarily induced by HCRV infection ( Figure 2A ). Remarkably, the differences in the transcriptome of N. benthamiana gradually increased with the time of HCRV infection ( Figure 2A ). The region of samples at 14 dpi of HCRV was separated the most from a region that contained the control samples, indicating that HCRV induced dramatic variations in the transcriptome of N. benthamiana at 14 dpi ( Figure 2A ). Further unsupervised hierarchical clustering of samples from the HCRV treatments reached a consensus with the related PCA result, showing that the main differences were generated by HCRV infection at 14 dpi ( Figure 2A ). Moreover, data showed that all the samples under the same treatments closely clustered on one branch, and a comparison of the difference between T14 vs. CK was more significant than that in T1 vs. CK and T7 vs. CK comparisons ( Figure 2A ). Additionally, the primary differences for TSWV infection represented by PC1 explained 36.9% of the variation between all the samples ( Figure 2B ). We found a similar phenomenon from the PCA and unsupervised hierarchical clustering analysis, which supported that the differences in the transcriptome of N. benthamiana gradually increased with the time of virus infection, while the primary differences were generated by virus infection at 14 dpi ( Figure 2B ). In contrast, for TH complex infection, the primary differences were generated by infection with the TH complex infection at 1 dpi, and there were 2,711 upregulated genes and 2,963 downregulated genes in the TH1 vs. CK comparison ( Figures 1C , 2C ). Both PCA and unsupervised hierarchical clustering analysis supported the credibility of our results that all the biological replicates were clustered together and separated from the other treatments ( Figure 2 ). Thus, we propose that the TH complex can infect plants in a different manner, which will help both viruses escape plant defense responses.

Subsequently, we performed qPCR assays to validate the gene expression values obtained by RNA-seq. For this purpose, we analyzed the expression of 10 tobacco resistance genes that are triggered by individual virus infection and suppressed by the TH complex ( Figure 2D ). We found that both qPCR assay and RNA-seq results reached a consensus that genes involved in JA-, MAPK-, and transcription factor-mediated defense responses and pathogenesis-related proteins, including NbPR1, NbPR4, and NbPR5, were triggered to upregulate by individual HCRV or TSWV, while these genes were suppressed to downregulate by the TH complex in vivo during the infection process ( Figure 2D ).

Considering the main variations in the transcriptome at 14 dpi of TSWV, the transcriptome data of tobacco in the T14 vs. CK comparison was used as a focus to investigate the transcriptomic changes associated with TSWV infection ( Figure 3 ). A total of 4,806 upregulated genes and 4,777 downregulated genes were identified as related to TSWV infection ( Figures 1C , 3A ). We then utilized MapMan software to visualize the involvement of these genes in plant responses to biotic stress ( Figure 3B ). As shown in Figure 3B , TSWV infection altered the expression pattern of the genes relevant to plant responses to biotic stress. We noted that TSWV primarily triggered the upregulation of genes involved in MAPK, WRKY, proteolysis, pathogenesis-related (PR) proteins, and secondary metabolism ( Figure 3B ). Most of the genes relevant to MAPK, proteolysis and WRKY were highly expressed in tobacco following TSWV infection at 14 dpi ( Figure 3B ). These genes have been well-characterized in plant defense responses to microbial infection, particularly MAPK, the PR proteins, and secondary metabolism (Pieterse et al., 2009; Zhou and Zhang, 2020). Thus, at 14 dpi, plants could deploy these resistance responses to defend against TSWV infection. Additionally, genes relevant to ethylene synthesis were upregulated in tobacco at 14 dpi, indicating the involvement of ETH in tobacco defense responses against TSWV ( Figure 3B ).

Subsequently, the DEs were collected for KEGG pathway and GO enrichment analyses ( Figures 3C, D ). The significant terms relevant to biological processes, molecular functions, and cellular components from the GO enrichment analysis are shown in Figure 3C . Among the terms relevant to biological processes, 306 terms were significantly enriched (FDR <0.05; Table S1 ), including photosynthesis, generation of precursor metabolites and energy, organic acid metabolic process, oxoacid metabolic process, carboxylic acid metabolic process, and alpha-amino acid metabolic process ( Figure 3C ; Table S1 ). The most significant terms for cellular components were plastid, chloroplast, thylakoid, chloroplast stroma, photosynthetic membrane, and plastid stroma, suggesting that TSWV infection primarily affected genes that function in plant photosynthesis and nutrient metabolism ( Figure 3C ; Table S1 ). Moreover, most of the genes that were categorized into molecular function terms were involved in cofactor binding, catalytic activity, isomerase activity, coenzyme binding, oxidoreductase activity, and kinase activity ( Figure 3C ; Table S1 ). Additionally, the KEGG pathway enrichment results showed that genes in several pathways were induced, including photosynthesis, biosynthesis of secondary metabolites, alpha-linolenic acid metabolism, and glycerolipid metabolism ( Figure 3C ). Over 80% of all the pathways shown in Figure 3D are involved in the plant metabolism of various essential metabolites, such as carbon metabolism, biosynthesis of amino acids, pyruvate metabolism, and glycerolipid metabolism ( Figure 3C ). Alterations in carbohydrate metabolism during the infection of pathogens that enable plants to reduce large amounts of starch granules in infected tissues to mobilize these storage molecules as complementary sources of energy have been well documented (Fan et al., 2010). Furthermore, resistance responses were significantly overrepresented in the DEGs, including secondary metabolism and the α-linoleic acid metabolism pathway associated with JA biosynthesis; most of the related DEGs were upregulated under TSWV infection ( Figures 3B, D ). Thus, we propose that reprogramming of these pathways is involved in the plant–TSWV interaction stage.

For HCRV infection, 3,272 upregulated genes and 3,253 downregulated genes from the H14 vs. CK comparison that responded to HCRV infection were functionally categorized according to Gene Ontology (GO) and KEGG pathway enrichment analyses ( Figure 4 ). The GO annotations of DEGs in tobacco following HCRV infection at 14 dpi were enriched for 10,267 GO terms, including 6,622, 2,510, and 1,135 terms relevant to biological process, molecular function, and cellular component, respectively ( Figure 4C ). In the “biological process” domain, the most frequent GO terms were photosynthesis, cofactor metabolic process, small molecule biosynthetic process, vitamin biosynthetic process, carboxylic acid biosynthetic process, carboxylic acid metabolic process, and oxoacid metabolic process ( Figure 4C ; Table S2 ). The results of the cellular component showed that most of the DEGs functioned in the plastid, chloroplast, thylakoid part, and plastid envelope ( Figure 4C ; Table S2 ), suggesting that HCRV also affected plant photosynthesis and essential metabolism. The molecular function GO category indicated that DEGs involved in lyase activity, transferase activity, carbon-carbon lyase activity, cofactor binding, catalytic activity, and kinase activity were highly represented ( Figure 4C ; Table S2 ). Consistent with the results observed with TSWV, the KEGG pathway enrichment results also showed the overrepresentation of pathways relevant to essential metabolite biosynthesis, including carbon metabolism, photosynthesis, biosynthesis of secondary metabolites, biosynthesis of amino acids, fructose and mannose metabolism, autophagy, and starch and sucrose metabolism ( Figure 4D ). In contrast, the significantly enriched pathway relevant to resistance responses was autophagy ( Figure 4D ). Autophagy has been reported to be involved in plant defense responses to eliminate viral infections (Chiramel et al., 2013). Although HCRV did not affect the biosynthesis of JA, it triggered the activation of secondary metabolism in plants, implying its conserved function in the responses of plants to viral infection ( Figure 4D ). Additionally, HCRV also induced changes in the expression of genes relevant to ethylene synthesis, MAPK, WRKY, PR proteins, and proteolysis ( Figure 4B ). Overall, we found that the responses of plants to TSWV and HCRV infections were highly conserved, and the TSWV and HCRV infections were closely associated with the processes of JA synthesis and autophagy, respectively.

Because the main changes between the TH vs. T and TH vs. H comparisons were at 7 dpi of infection, we further analyzed the transcriptome profiles of tobacco at 7 dpi of HCRV, TSWV, and TH complex infections. In this study, we identified 6,227 upregulated genes and 2,445 downregulated genes in the H7 vs. TH7 comparison, whereas the T7 vs. TH7 comparison had 6,575 upregulated genes and 4,013 downregulated genes ( Figure 5A ). We then analyzed these DEGs from both pairwise comparisons by KEGG pathway and GO enrichment analyses. A functional analysis showed that both sets of DEGs from the TH vs. H and TH vs. T comparisons were enriched in 9,542 and 10,133 GO terms, respectively ( Tables S3 , S4 ). The results of GO enrichment on the DEGs from the H7 vs. TH7 comparison showed that these genes were highly enriched for several processes, including rRNA binding, lyase activity, peptide biosynthetic process, photosynthesis, carboxylic acid metabolic process, plastid, chloroplast, and ribosome ( Table S3 ). Moreover, both enriched GO results from TH7 vs. T7 highly overlapped with those from the TH14 vs. H14 comparison ( Tables S3 , S4 ). GO enrichment analysis on the DEGs from the T7 vs. TH7 comparison displayed the overrepresentation of plastid, chloroplast, photosynthesis, carboxylic acid metabolic process, catalytic activity, isomerase activity, and lyase activity ( Table S3 ). Both results suggested that complex infections with HCRV and TSWV could primarily alter plant photosynthesis and metabolism for essential metabolites and nutrients.

Subsequently, a pathway analysis was performed to predict the significantly enriched pathways that involved the DEGs from TH7 vs. T7 and TH7 vs. H7 comparisons to investigate the differences between individual infection of HCRV or TSWV and the TH complex infection. The significantly enriched pathways in the TH7 vs. T7 comparison were photosynthesis, ribosome, carbon metabolism, biosynthesis of secondary metabolites, biosynthesis of amino acids, starch and sucrose metabolism, and plant hormone signal transduction ( Figure 5C ). Moreover, the data from the TH7 vs. H7 comparison showed that the differences were primarily generated in the ribosome, photosynthesis, biosynthesis of amino acids, carbon metabolism, and biosynthesis of secondary metabolites ( Figure 5B ). Numerous pathways of essential metabolism were affected by TH infection at 7 dpi compared with neither HCRV nor TSWV ( Figures 5B, C ). Remarkably, some of these pathways are known to play an important role in plant defense responses to HCRV or TSWV infections, particularly plant hormone signal transduction and biosynthesis of secondary metabolites ( Figures 5B, C ). Importantly, we noted that many of the genes involved in plant basal resistance, particularly JA biosynthesis, MAPK, secondary metabolism, and pathogenesis-related proteins, could not be activated by TH complex infection but were triggered to upregulate under HCRV and TSWV infections ( Figure 5D ). Among pathogenesis-related proteins, we identified significant downregulation of NbPR1, NbPR4, and NbPR5 in tobacco under TH complex infection compared to individual virus infection ( Figure 2D ). Then, we constructed tobacco with the silencing of these PRs to investigate their involvement in tobacco resistance to HCRV or TSWV. RT-qPCR assay confirmed the silencing of these genes in representative tobacco mutants ( Figure 5E ). Typically, NbPR1 expression was suppressed by about 79% in tobacco compared to the controls ( Figure 5E ). Then, HCRV or TSWV were inoculated on these tobaccos. Fourteen days post-inoculation, we found observable symptoms caused by HCRV or TSWV on tobacco leaves compared to the controls ( Figure 5E ). Then, we determined the biomass of HCRV or TSWV in related tobacco, and the results showed that the biomass of HCRV or TSWV in tobacco with NbPR1 silencing was 6.4- and 4.7-fold higher than that in the controls ( Figure 5E ), suggesting that silencing of NbPR1 resulted in loss of resistance in tobacco. Thus, we propose that the TH complex can synergistically infect plant tissue to escape or suppress the plant’s basal resistance for their successful infection.

MicroRNAs have been found to play important roles in plant–virus interactions, and they could interact with their target genes to degrade related mRNA or suppress its translation to regulate the series of plant responses (Ding, 2010; Finnegan and Pasquinelli, 2013; Borges and Martienssen, 2015). In this study, we further detected the dynamic expression variations of the miRNAs of tobacco following HT complex infection at 1, 7, and 14 dpi ( Figure 6A ). A total of 2,952 miRNAs were identified from all the sample libraries. In total, 517, 132, and 222 miRNAs were differentially expressed in the TH1 vs. CK, TH7 vs. CK, and TH14 vs. CK comparisons, respectively ( Figure 6A ). In addition, we found that a total of 346 miRNAs were upregulated in tobacco at 1, 7, and 14 dpi of TH infection, suggesting their involvement in plant–TH communication during the development of TH infection ( Figures 6A, D ). We then constructed the co-expression network of these 346 miRNAs with all the mRNAs of tobacco ( Figure 6B ). In addition, all the negative relationships between miRNA and mRNA were used to construct the topological network as shown in Figure 6B . Remarkably, 28 miRNAs were identified as hub regulators in this topological network, including 19 upregulated miRNAs at 1 dpi of TH infection and nine miRNAs at 7 dpi of TH infection ( Figure 6B ; Table S4 ). We then collected all the genes that negatively correlated with these miRNAs in the network and performed KEGG pathway enrichment analysis on these genes ( Figure 6C ). The results showed that these genes are primarily involved in plant–pathogen interaction, terpenoid biosynthesis, phagosomes, MAPK signaling, alpha-linolenic acid metabolism, and plant hormone signal transduction ( Figure 6C ). The most significant pathways were plant–pathogen interaction, MAPK signaling, linoleic acid metabolism, and metabolisms for secondary metabolites ( Figure 6C ), and these pathways that were significantly enriched were resistance-related processes associated with HCRV and TSWV infection, suggesting that the complex could trigger miRNAs to suppress plant basal resistance. Among these miRNAs, we noted that they exhibited two patterns of expression. A total of 19 miRNAs were highly expressed in tobacco at 1 dpi of TH complex infection, while nine miRNAs were primarily upregulated in TH-challenged tobacco at 7 dpi ( Figure 6D ).

To investigate the effect of these miRNAs associated with TH complex infection on mediating plant resistance, we overexpressed the top five miRNAs from the network, including novel-m0693-5p, novel-m0782-3p, novel-m0915-3p, novel-m1488-5p, and novel-m1992-3p in tobacco. Consistently, RT-qPCR assay confirmed that the TH complex triggered the upregulation of these five miRNAs, while individual HCRV or TSWV did not affect their expression or even suppress their expression in tobacco during the infection process ( Figure 7A ). Then, the expression levels of resistance genes that were triggered by infection with individual HCRV or TSWV viruses were further detected using RT-qPCR. The results showed that these resistance genes relevant to JA-, calcium-, and WRKY-mediated defense were dramatically inhibited by overexpression of these miRNAs ( Figure 7B ). Typically, all these miRNAs significantly inhibited the expression of LOX6, while novel-m0782-3p and novel-m0915-3p inhibited AOC expression in tobacco ( Figure 7B ), which could further suppress JA biosynthesis and related defense responses to reduce resistance in tobacco to virus. For calcium-related genes, individual novel-m0693-5p, novel-m0915-3p, or novel-m1992-3p overexpression mainly inhibited the expression of CPK18 and CPK4 that are involved in regulating calcium-related defense responses ( Figure 7B ). Additionally, all these miRNAs inhibited the expression of WRKY6, while novel-m0693-5p, novel-m0915-3p, or novel-m1992-3p overexpression also inhibited WRKY1 and WRKY22 in vivo ( Figure 7B ). Additionally, miR1992-3p significantly inhibited the expression of NbSAM in vivo ( Figure 7B ). Taken together, these results suggest the involvement of these 28 miRNAs, especially novel-m0693-5p, novel-m0782-3p, novel-m0915-3p, novel-m1488-5p, and novel-m1992-3p, in the early communication between the plant and virus, which could further suppress the host defense responses so that both HCRV and TSWV could successfully infect the plant.

It is well known that miRNAs exert their functions by inhibiting the expression of target genes, and miRNAs need to be strictly complementary with their target genes to cleave at the pairing sites (Zhao et al., 2015). To reveal the target of these miRNAs, we predicted their potential targets through the website psRNATarget (http://plantgrn.noble.org/psRNATarget/). Under relatively strict parameters, NbWRKY6 and NBSAM were predicted as the putative targets of novel-m0782-3p and miR1992-3p, respectively. Of these predicted target genes, we used 5’ RLM-RACE to validate the miRNA cleavage sites in target mRNAs. We found cleavage sites for novel-m0782-3p and miR1992-3p in NbWRKY6 and NBSAM ( Figures 7C, D ), which supports the view that these genes are the direct targets of the corresponding miRNAs in tobacco. Among these target genes, a cleavage site for novel-m0782-3p was detected in the predicted miRNA complementary region, while for miR1992-3p, cleavage sites were detected in the predicted miRNA complementary region and upstream of this region ( Figure 7B ). Although NbPR1 was not directly cleaved by novel-m0782-3p and miR1992-3p, we found that silencing of NbWRKY and NbSAM could result in downregulation of NbPR1 ( Figure 7E ), suggesting that NbWRKY and NbSAM may be involved in regulating NbPR1. Overall, these results suggested that the TH complex could trigger novel-m0782-3p and miR1992-3p to suppress NbWRKY6 and NBSAM expression, leading to downregulation of NbPR1 and decreases in tobacco resistance to both viruses.