Leaves and twigs of gardenia showing mosaic, mottle, or chlorotic symptoms were collected from a small garden at Ningbo University in Ningbo City, Zhejiang Province, China. After collection, samples were rapidly frozen in liquid nitrogen to preserve RNA integrity. Following grinding, approximately 50 mg of the ground powder underwent total RNA extraction with the TRIzol reagent (Invitrogen, Carlsbad, CA, USA). After the assessment of RNA quality and quantity, RNAprep Pure Plant Kit (Polysaccharides & Polyphenolics-rich) (TianGene, Beijing, China) was employed for total RNA extraction, and the samples were subsequently sent to Novogene (Beijing, China) for both metatranscriptome (RNA sequencing [RNA-Seq]) and small RNAs sequencing.

The sequencing library preparation and sequencing were carried out by Novogene. In brief, for RNA-Seq, the RNA sequencing library was prepared by deleting ribosomal RNAs from the total RNAs using the Ribo-Zero rRNA Removal Kit (Plant) (Illumina, San Diego, CA, USA), following the manufacturer’s instructions. Subsequently, the Illumina Novaseq 6000 system was employed for sequencing in a paired-end mode with a read length of 150 bp on both ends (150 bp × 2). For small RNA sequencing, the small RNA sequencing library was prepared with TrueSeq Small RNA Library Prep Kit (Illumina) and the sequencing was carried out on an Illumina HiSeq 4000 system with a read length of 50 bp.

Raw RNA-Seq dataset was processed to remove adapter regions and low-quality reads using fastp (Chen et al., 2018). Reads shorter than 36 nucleotides were removed. De novo assembly of the remaining reads were conducted using Trinity assembler (version 2.8.5) (Hass et al., 2013). To identify and annotate the virus-associated contigs, the assembled transcriptome contigs were subjected to BLASTx/BLASTn analysis against the NCBI (https://www.ncbi.nlm.nih.gov) NR and NT databases (accessed on May 2023) (Altschul et al., 1990). Additionally, a prebuilt RdRP Hidden Markov Model database was searched using the hmmsearch tool embedded in HMMER (http://hmmer.org/) to detect highly divergent viruses which may exist.

Total RNA was extracted from leaf samples of eight plants, including five symptomatic plants and three asymptomatic plants, using TRIzol reagent (Invitrogen) according to the manufacturer’s instructions. To confirm the infection of potyvirus, the potyvirus degenerate primers PotyF and PotyR (Marie-Jeanne et al., 2000) were used in RT-PCR. Subsequently, to recover the full length of BCMV, the whole genome regions were amplified by overlapping RT-PCR and RACE. Each amplicon was inserted into the pClone007 Blunt Simple Vector (Tsingke, Beijing, China) and sequenced by Sanger sequencing. The primers used in this study are listed in Supplementary Table 1.

To confirm the evolutionary status of BCMV-gardenia in the genus Potyvirus, genome sequences of representative members within the genus Potyvirus were retrieved from NCBI GenBank database. The polyprotein encoding regions of these viruses were extracted and aligned using MAFFT (Katoh and Standley, 2013). A phylogenetic tree deduced based on the maximum likelihood method was constructed using IQ-TREE (v1.6.6) (Nguyen et al., 2015). The nucleotide substitution model GTR + F + I + G4, selected by ModelFinder (Kalyaanamoorthy et al., 2017), was applied. Confidence in the topology was assessed using 1,000 bootstrap replicates.

To study the evolutionary relationship among different isolates of BCMV identified in various hosts, the BCMV isolate sequences were retrieved from NCBI GenBank database and similar evolutionary analysis was carried out with the best-fitted model selected by ModelFinder.

The adapter sequence was removed using the cutadapt tool (v1.16) (Martin, 2011). Reads that were 18 nt to 30 nt and contained no ambiguous nucleotides (‘N’) were retained as clean reads. These clean reads were then clustered into unique reads using an in-house Python script. The remaining small RNA clean reads were subsequently aligned to the BCMV-gardenia virus reference genome (accession no. OR832357) using the bowtie tool (Langmead et al., 2009), allowing for up to one mismatch. The profile of virus-derived small interfering RNA (vsiRNAs), including length distribution, 5′ terminal nucleotide preference, genome distribution, and polarity, was explored through the analysis of small RNAs that could be mapped to the virus reference genome.

To elucidate the potential pathogens of gardenia exhibiting disease symptoms, leaves and twigs showing symptoms such as mosaic, chlorosis were collected from a small garden in Ningbo University, Ningbo City, Zhejiang Province, China (Fig. 1). About 5-7 different leaves were selected randomly to extract the total RNA for both RNA-Seq and small RNA sequencing. After RNA-Seq, a total of 43,103,034 reads were produced. Following quality control, more than 41,319,142 reads (95.86%) were retained for sequence assembly. For small RNA sequencing dataset, after strict quality control, only 14,807,281 of the 30,234,844 raw reads remained for subsequent analysis (Table 1).

Then the clean reads in RNA-Seq library were assembled to long contigs and annotated to NCBI NR/NT reference database. A large contig with a length of 10,043 nt was annotated to a BCMV isolate genome sequence (KF114860) with significant similarity (91.66% nucleotide similarity in BLASTn). A large polyprotein sequence predicted from that viral contig, 3,222 aa in length, could be annotated to polyprotein of BCMV (AYE54419) in NCBI-BLASTp with 96.61% similarity. This large contig also could be annotated to contain a potyviridae domain with e-value of 5.99e-230 in hmmsearch analysis. Moreover, both BLASTx search-based and hmmsearch-based methods showed that there were no other known or unknown viruses in the RNA-Seq library. According to species demarcation criteria of the family Potyviridae proposed by International Committee on Taxonomy of Viruses (<76% nucleotide identity and <82% amino acid identity in complete ORF sequence), the viral contig we identified was considered as a strain of BCMV and named BCMV-gardenia.

To confirm the presence of the potyvirus, the potyvirus degenerate primers PotyF and PotyR (Marie-Jeanne et al., 2000) were synthesized and used in RT-PCR. After PCR amplification, the bright bands with expected length (~300 bp) were only detected in samples from five symptomatic leaves and not in three asymptomatic leaves (Fig. 2). The sequences of amplicons showed approximately 95.41% nucleotide sequence identity to BCMV isolate Ir-MSC1 (MF498888). Subsequently, the complete nucleotide sequences of this BCMV candidate virus were amplified by overlapping RT-PCR and RACE. The gel electrophoresis image depicted in Supplementary Fig. 1 displayed that all these sets of specific primers produced specific bands with corresponding length (Supplememtary Table 1).

After Sanger sequencing of the amplification product from the overlapping RT-PCR and RACE, the full-length genome sequence of BCMV-gardenia was found to be 10,054 nt in length, excluding the poly(A) at the 3′ termini, with a 5′ untranslated region (UTR) of 134 nt and a 3′ UTR of 286 nt (GenBank accession no. OR832357) (Fig. 3A). Consistent with other potyviruses, nine highly conserved proteolytic cleavage sites were predicted, and ten putative mature proteins could be identified, including P1 (443 aa, 49.96 kDa), HC-Pro (457 aa, 52.06 kDa), P3 (347 aa, 39.57 kDa), 6K1 (52 aa, 5.93 kDa), CI (634 aa, 71.25 kDa), 6K2 (53 aa, 6.15 kDa), VPg (190 aa, 21.90 kDa), NIa-Pro (243 aa, 27.56 kDa), NIb (516 aa, 59.62 kDa), and CP (287 aa, 32.20 kDa) from N-termini to C-termini. The putative P3N-PIPO, encoding a 74 aa protein with 8.77 kDa, was also identified at nucleotide position from nucleotide position 3,292 nt to 3,514 nt (Fig. 3A and C). To compare both the similarities and variation in the conserved protein cleavage sites among different BCMV isolates, as well as bean common mosaic necrosis virus (BCMNV) and SMV, all belonging to the genus Potyvirus, the conserved protein cleavage sites were predicted and indicated in the Fig. 3A. From the results, it was evident that despite differences in the polyprotein sequence between BCMV-gardenia and BCMV with accession number NC003397 (identified in black-eye cowpea), their conserved protein cleavage sites exhibit complete identity. For BCMV isolates, BCMNV, and SMV, differences were only observed in the protein cleavage sites between P1 and HC-Pro proteins (Fig. 3A). We also compared the amino acid and nucleotide similarity of each protein between BCMV-gardenia and another BCMV isolate (NC003397). We found that these two BCMV isolates showed higher similarity at the C-termini of the polyprotein, such as CP protein region, while the sequence similarity is lowest at the N-termini, especially in the P1 protein region (Fig. 3C). When clean reads were mapped back to the BCMV-gardenia reference genome, we found that every region of the viral genome was covered, with the highest depth at both ends, especially the 3′ end of the viral genome (Fig. 3B).

To confirm the evolutionary status of BCMV-gardenia in the genus Potyvirus, the phylogenetic tree of BCMV-gardenia and other representative member of the genus Potyvirus, was constructed by a maximum likelihood algorithm using IQ-TREE program. The multiple sequence alignment was performed by MAFFT. As shown in Supplementary Fig. 2, BCMV-gardenia was highly associated and clustered with BCMV isolate (NC003397), reinforcing the conclusion that virus identified in this research is an isolate of BCMV. To study the evolutionary relationships among different BCMV isolates, a midpoint-rooted phylogenetic tree was constructed using the various isolate of BCMV identified in different hosts (Fig. 4). Except for the BCMV-gardenia, the genome of other BCMV isolates were retrieved from NCBI GenBank database. The final dataset used to construct evolutionary tree comprised 101 sequences, the vast majority of which represented the BCMV genomes identified in leguminous plant hosts. However, there were three exceptions: BCMV-gardenia (identified in gardenia), KM076650 (identified in Cudrania tricuspidata), and MZ670770 (identified in Nandina domestica). From the evolutionary tree, we found that BCMV isolate identified in the same host species tend to cluster into a phylogenetic group. Remarkably, the three BCMV isolates identified outside leguminous hosts, were adjacent to each other on the whole evolutionary tree (Fig. 4). Notably, among the three BCMV isolate identified outside the leguminous hosts, KM076650 (identified in Cudrania tricuspidata) and MZ670770 (identified in Nandina domestica) formed a monophyletic group and BCMV-gardenia (identified in gardenia) seemed to have a closer evolutionary relationship with other BCMV identified in leguminous plants (Fig. 4).

RNA silencing serves as a highly effective defense mechanism against viral infections employed by plants. This defense strategy relies on the production of the vsiRNAs. Therefore, the characterization of the vsiRNAs profile contributes to our understanding of the virus-host interaction. We found that the vsiRNAs derived from BCMV-gardenia are mainly 21nt and 22nt, with 21nt being the most abundant (Fig. 5A). We investigated the 5′-terminal nucleotide preference of vsiRNAs derived from BCMV-gardenia. Our analysis revealed a predominant U/A preference for vsiRNAs with lengths of 21 nt and 22 nt (Fig. 5A) or across all vsiRNAs (Fig. 5B). Regarding the genomic regions from which vsiRNAs were generated, our observations indicated a relatively even distribution across the entire viral genome. However, specific local regions emerged as hot spots for vsiRNA production, with 21 nt and 22 nt vsiRNAs being the predominant types in each region (Fig. 5C). We further explored the strand polarity and specific origin sites of vsiRNAs along BCMV-gardenia in detail. The results indicated that 51.1% of vsiRNAs originated from the sense strand, while 48.9% originated from the antisense strand (Fig. 5D).