Hibiscus rosasinensis were sampled from Taiwan Farmers Entrepreneurship Park in Zhangzhou, Fujian. Total RNA was extracted from both healthy and virus-infected leaves of Hibiscus rosa-sinensis with Trizol reagent (Invitrogen, Carlsbad, CA, USA) as instruction. Concentrations and integrity of total RNA were detected using a spectrophotometer (Nanodrop 2000, Thermo Fisher Scientific, Wilmington, DE, USA) and a BioAnalyser 2100 (Agilent, Palo Alto, CA, USA), respectively. To detect infection of HLFPV, primers (available upon requested) were designed and synthesized in SANGON BIOTECH Co., Ltd. (Shanghai, Beijing, China) basing on the sequences of coat protein (CP) gene of the HLFPV Japanese isolate (Yoshida et al., 2014). Reverse transcription polymerase chain reaction (RT-PCR) was conducted with the primers to amplify the CP gene with FastKing One Step RT-PCR Kit (TIANGEN BIOTECH, Beijing, China) as instruction. To characterize the sequences properties of HLFPV genome, 11 primers pairs (available upon requested) for amplifying 11 sequence segments with overlap regions were designed and synthesized in SANGON BIOTECH Co., Ltd. basing on the HLFPV Japanese isolate (Yoshida et al., 2014). Then, 11 products of RT-PCR were cloned into pMD18-T vector. The obtained recombined vectors were sequenced and the obtained sequences were combined into a full sequence with software DNAMAN. ORFs were predicted on the NCBI website, https://www.ncbi.nlm.nih.gov/orffinder/. The phylogenetic tree was established by maximum likelihood method with MEGA 5.1 software with 1,000 replicates. Sequence similarity comparison of HLFPV ORFs was performed with NCBI Blast. Sequence similarity comparison of HLFPV 5′-terminal and 3′-terminal regions was performed with software BioEdit. Prediction of secondary structures of the full-length viral genome and its 3′-terminal UTR was performed with RNAfold server at the website, http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi.

Leaf dips of both healthy and virusinfected leaves of hibiscus plants, Hibiscus rosa-sinensis were made on 200 mesh Formvarcoated copper grids and stained with 5% (w/v) uranyl acetate. Virus particles were observed under H-7650 Hitachi transmission electron microscope (Hitachi, Tokyo, Japan) at 80 kV as described previously (Lan et al., 2016).

Total RNA was extracted from virus-infected and healthy leaves of hibiscus plants, Hibiscus rosa-sinensis using Trizol Reagent (Invitrogen) as introduction. Concentrations and integrity of total RNA were detected using a spectrophotometer (Nanodrop 2000) and a BioAnalyser 2100 (Agilent Technologies), respectively. Small RNA libraries were constructed and sequenced with the Solexa protocol (Gen Denovo, Guangzhou, China) as described previously (Lan et al., 2015, 2016). Briefly, small RNA molecules (18-32 nt) were purified and Solexa adaptors were ligated to 5′-and 3′-terminals. Then, these small RNA molecules were amplified with the adaptor primers for 25 cycles and fragments of about 90 bp (small RNA and adaptors) were isolated from the agarose gel. The purified DNA was utilized directly for small RNA sequencing analysis with Illumina’s Solexa Sequencer. Raw data sets for the small RNA were analyzed. In brief, adaptor sequences were trimmed and small RNA reads without an identifiable linker were removed. The remaining reads were filtered by length. Reads of >32 or <18 nt were discarded. To identify vsiRNAs, using the software Bowtie v.12.7 with a parameter of 0 mismatch, we aligned all the cleaned reads to the HLFPV genome. The downstream analyses of vsiRNAs were performed using Perl scripts and Excel. Target genes were predicted with the small RNA target analysis server Patmatch_v1.2 with a parameter of 0 mismatch, minimum free energy ≥90% of free energy when vsiRNAs perfectly binds to target genes. The Gene Ontology (GO) annotation analysis of the predicted target genes of the vsiRNAs was performed using BLAST2GO software (Lan et al., 2019).

Hibiscus plants samples, Hibiscus rosa-sinensis showing severe chlorotic spots, ringspots and chlorotic mottle on plant leaves (Fig. 1A) for about 10 days were collected from Taiwan Farmers Entrepreneurship Park in Zhangzhou, Fujian. The infection of HLFPV was concurrently confirmed by RT-PCR and by H-7650 Hitachi transmission electron microscopy. Our results showed that CP genes (about 500 bp) of HLFPV were successfully amplified by RT-PCR with specific primers (available upon requested) from plant leaves showing symptom but not healthy plant leaves (Fig. 1B). The obtained positive RT-PCR products were sequenced and the results of NCBI Nucleotide Blast confirmed that it matched with CP gene of HLFPV genome (data not shown). The phylogenetic tree of the CP gene showed that HLFPV-Ch isolate clustered with analogue of other HLFPV isolates (Supplementary Fig. 1), implying high homology among them. Furthermore, electron microscopy by negative staining method revealed that approximately 300 nm rod-shaped virus particles, perhaps representing for HLFPV, existed in leaf dips prepared from viruses-infected leaves of Hibiscus rosa-sinensis (Fig. 1C) but not from healthy leaves of Hibiscus rosa-sinensis (data not shown). Thus, these results preliminary proved Hibiscus rosa-sinensis may be infected by HLFPV. This is reported for HLFPV infection in China for the first time.

By being cloned, sequenced and joined, the complete genomic sequence of HLFPV isolated from China was 6,465 nucleotides (nt) in length (Fig. 2A). It contained four ORFs encoding proteins of 128.5, 186.5, 31.1 and 17.6 kDa in ORF1 (nt 47-3,463), ORF2 (nt 47-4,987), ORF3 (nt 4,968-5,822), and ORF4 (nt 5,794-6,270), respectively (Fig. 2A) by being predicted on the NCBI website. ORF3 overlaps partially with ORF2 and ORF4 by 20 and 29 nt, respectively (Fig. 2A). These results implied the similarity in genomic organization or arrangement of tobamoviruses. The four ORFs of HLFPV isolated from China were 99.4-100%, 99.4-100%, 99.3-100%, and 99.2-100% identical at the nucleotide level and 99.4-100%, 99.4-100%, 99.7-100%, and 98.7-100% identical at the amino acid level to counterparts of other four HLFPV isolates, respectively (Table 1). These results suggested that the detected HLFPV isolated from China in Hibiscus rosa-sinensis was a tobamovirus. Subsequently, we designated this China isolate of HLFPV as HLFPV-Ch.

The length of the 5′-ternimal UTR of HLFPV-Ch was 46 nt, which is shorter than HLFPV Japanese and Brazil isolates but longer than Italy isolate (Gao et al., 2016; Nerva et al., 2018; Yoshida et al., 2014). The 5′-UTR shared 21.7-43.7% identity at the nucleotide level with other HLFPV isolates and showed the highest identity to that of HLFPV Italy isolate (Table 1). Eight copies of CAA triplets, which are conserved in the 5′-ternimal UTR regions of other tobamoviruses, were found in the 5′-ternimal UTR of HLFPV-Ch. Furthermore, the 5′-proximal 10 nucleotides of all HLFPV isolates including HLFPV-Ch was GUAUGUUUUA, implying its importance in the infection cycles of virus.

The 3′-ternimal UTR of HLFPV-Ch was 195 nt in length. It was 36.9-49.7% identical at the nucleotide level to counterparts of other three HLFPV isolates, showing the highest (49.7%) identity to that of HLFPV isolated from Brazil (Table 1). The 3′-ternimal UTR of HLFPV-Ch contained a 79 nt poly(A) region but not the pseudoknot region (Fig. 2A), both of which was conserved in most tobammoviruses (Yoshida et al., 2014), at the downstream of the stop codon in ORF4. More and more evidences have shown that poly(A) region of plant viruses plays an important role in the infection cycles of viruses by interactions with many host proteins, such as poly(A)-binding proteins (Beauchemin and Laliberté, 2007). It has been shown that downstream of the poly(A) region of the 3′-UTR of HLFPV isolated from Japan was predicted to contain a tRNA-like secondary structure (Yoshida et al., 2014). Using RNAfold web server (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi), downstream of the poly(A) region of the 3′-UTR of HLFPV-Ch was also predicted to contain a tRNA-like secondary structure (Supplementary Fig. 2), implying their the similar functions in virus infection cycle. The 3′-proximal sequence of HLFPV-Ch was also CCCA, suggesting the sequence conservatism among tobamoviruses.

Following determination of HLFPV-Ch infection, small RNA libraries were constructed and sequenced with the Solexa protocol (Gen Denovo) with total RNA extracted from HLFPV-Ch-infected and healthy Hibiscus rosa-sinensis. A total of 10,839,244 and 8,733,689 sRNAs reads were sequenced from HLFPV-Ch-infected and healthy Hibiscus rosa-sinensis, respectively. It showed that 21-, 24-, and 22-nt sRNAs were most dominant in HLFPV-Ch-infected and healthy Hibiscus rosasinensis (Fig. 3A). Furthermore, both 21-nt, 22-nt and 23-nt sRNAs reads increased markedly, but 24-nt sRNA reads decreased significantly in HLFPV-Ch-infected Hibiscus rosa-sinensis (Fig. 3A), which is similar to previous reports for the same Tobamovirus genus (Lan et al., 2019; Li et al., 2016). These results suggested that the infection of virus, HLFPV-Ch, modulated the distribution pattern of sRNAs in Hibiscus rosa-sinensis.

After mapped to the HLFPV-Ch genome, 30,925 unique vsiRNAs (18-32 nt) were identified in HLFPV-Ch-infected Hibiscus rosa-sinensis, accounting for 0.29% of siRNAs reads. Conversely, only 571 unique vsiRNAs (18-32 nt), accounting for 0.007% of siRNAs reads, were identified in healthy Hibiscus rosa-sinensis. In HLFPV-Ch-infected Hibiscus rosa-sinensis, the majority (71%) of vsiRNAs were 19-24 nt in length, with relative average abundance for different sized species although 21-nt and 22-nt vsiRNAs being most abundant orderly (Fig. 3B); however, the majority of vsiRNAs in many virus-plant systems were 21-nt and 22-nt in length (Kreuze et al., 2009; Lan et al., 2019; Li et al., 2016; Mitter et al., 2013; Xu and Zhou, 2017; Yan et al., 2010; Yang et al., 2014). Thus, these results implicated that the homologues of DCL4 and DCL2 in Hibiscus rosa-sinensis may also be the predominant DCLs involved in the biogenesis of 21-nt and 22-nt vsiRNAs respectively which functions as the predominant antiviral silencing mediators (Blevins et al., 2006; Deleris et al., 2006; Ding, 2010; Donaire et al., 2008; Liu et al., 2018; Niu et al., 2017); as for other DCLs in Hibiscus rosa-sinensis participating in different-lengths vsiRNAs generation will be the future discussed topic. Furthermore, the distance of 5′ ends (overhangs) of 21-nt vsiRNAs duplexes were analyzed. The results suggested that 21-nt vsiRNA duplexes with 2-nt overhangs were the most abundant 21-nt duplexes, followed by 1-nt overhangs and then 0-nt overhangs 21-nt duplexes in HLFPV-Ch-infected Hibiscus rosa-sinensis (Fig. 3C), which was similar that of other viruses-plant and viruses-invertebrates systems (Liu et al., 2018; Niu et al., 2017). Thus, the results revealed that 21-nt vsiRNAs duplexes with 2-nt overhangs were the most efficient triggers of RNAi in Hibiscus rosa-sinensis.

Many researches have suggested that the 5′-terminal nucleotides of vsiRNAs have a important role in modulating the sorting of vsiRNAs to different AGO complexes in plants (Mi et al., 2008). In this study, our analysis showed that HLFPV-Ch unique vsiRNAs demonstrated a clear tendency to begin with uracil (U, 31%) and adenine (A, 28%) as compared with cytosine (C, 23%) and guanidine (G, 17%) (Fig. 3D), which was consistent with previous studies for diverse plant-virus systems (Donaire et al., 2008, 2009; Lan et al., 2019; Xu and Zhou, 2017). To obtain further understanding of HLFPV-Ch vsiRNA sorting, the complexity of vsiRNAs was analyzed for different sizes species. For 21-nt and 22-nt HLFPV-Ch vsiRNAs, a clear preference for vsiRNAs beginning with a 5′-U was observed (Fig. 3D). This is consistent with the function described for AGO1 in defending against RNA viruses in plant (Morel et al., 2002; Qu et al., 2008). However, for 23-nt, 24-nt, and 25-nt HLFPV-Ch vsiRNAs, a strong bias for vsiRNAs beginning with a 5′-A was observed (Fig. 3D), implying the high binding affinity of AGO2 and AGO4 for these vsiRNAs. The low proportion of vsiRNAs beginning with G was also observed and consistent with previous reports of pathogen-systems (Donaire et al., 2009; Li et al., 2016; Mi et al., 2008; Mitter et al., 2013). Thus, our finding suggested the involvement of different AGO proteins in binding of vsiRNAs in different plant species. In future, isolation and determination of AGO-vsiRNAs complexes from hibiscus plants, Hibiscus rosa-sinensis infected with HLFPV-Ch will provide us more functions information about vsiRNAs sorting in hibiscus plants.

To reveal the origin of vsiRNAs, we aligned vsiRNAs with positive and negative strands of HLFPV-Ch genome. Our results showed that HLFPV-Ch vsiRNAs were derived predominantly from the viral positive-strand RNA, accounting for about 86% (Fig. 3E). This is broadly contrasted to many ssRNA positive-strand viruses, including the cucumber green mottle mosaic virus of the same genus, vsiRNAs of which were produced nearly equally from the positive and the negative strands (Chen et al., 2018; Ho et al., 2007; Lan et al., 2020; Li et al., 2016; Xia et al., 2014; Yang et al., 2014). Additionally, single-base resolution maps of total unique vsiRNAs along with HLFPV-Ch genomes were created using Bowtie tools and in-house Perl scripts. Our results showed that the majority vsiRNAs were derived from the intragenic regions, not the 5′-or 3′-terminal regions of HLFPV-Ch genome (Fig. 2B). These vsiRNAs had a continuous but heterogeneous (Hot spot and Cold spot) distribution along the HLFPV-Ch genome (Fig. 2B). Taken together, our results implied that highly structured regions present in HLFPV-Ch genome may be the substrates for DCL enzymes to cutting into vsiRNAs. To confirm the speculation, we evaluated the secondary structures of HLFPV-Ch genome with RNAfold server. Our results showed that lots of stem loop secondary structures were formed in intragenic regions and corresponded to the hotspots (Fig. 2C and D). Taken together, these results supported that vsiRNAs should originate predominantly by direct DCL cleavage of imperfect duplexes in the most folded regions of the positive strand of HLFPV-Ch RNA molecular in Hibiscus rosa-sinensis. This result is similar to the speculation for vsiRNAs origin from plants infected with diverse ssRNA positive-strand viruses (Molnár et al., 2005).

Identification and characterization of target genes is very important to understand deeply the functions and implications of vsiRNAs in interaction with host. In this study, we used the small RNA target analysis server Patmatch_v1.2 to predict putative target genes of Hibiscus rosa-sinensis by HLFPV-Ch-derived vsiRNAs. Totally, 722 putative target genes were predicted by HLFPV-Ch-derived siRNAs (Supplementary Table 1). In BLAST2GO analysis, these predicted target genes are involved in abroad range of cellular component, molecular function and biological processes (Fig. 4, Supplementary Table 1). The three most highly represented GO terms were “intracellular membrane-bounded organelle”, “membrance” and “intracellular cytoplasm” under the cellular component category, comprising 90% of the annotated genes (Fig. 4A). The two most abundant GO terms under the molecular function category, “catalytic enzyme activity” and “binding”, represented 90% of the annotated genes (Fig. 4B). The terms “metabolic process”represented 73% of the classified targets under the biological process section (Fig. 4C). These results suggested that vsiRNAs may play vital roles due to their interactions with host by targeting their genes during HLFPV-Ch infection in hibiscus plants, Hibiscus rosa-sinensis.