H. rosa-sinensis showing leaf mottle and chlorotic spots symptoms were sampled from Taiwan Farmers Entrepreneurship Park in Zhangzhou city in Fujian province of China. Total RNA was extracted from both healthy and virus-infected leaves of H. 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, Waltham, MA, USA) and a BioAnalyser 2100 (Agilent Technologies, Waldbronn, Germany), respectively. To detect HCRSV infection, primers (available upon requested) were designed and synthesized in Sangon Biotech Co., Ltd. (Shanghai, China) based on the sequences of coat protein (CP) gene of the HCRSV Fujian isolate (Zheng et al., 2018). 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 HCRSV genome isolated from Zhangzhou city, Fujian province, China (named HCRSV-ZZ), seven primers pairs (available upon requested) for amplifying seven sequence segments with overlap regions were designed and synthesized in Sangon Biotech Co., Ltd. basing on the genome of HCRSV Fujian isolate (Zheng et al., 2018). Then, these seven products of RT-PCR were cloned into pMD18-T vector. The obtained recombined vectors were sequenced and the obtained sequences were combined into a relatively longer sequence with software DNAMAN. Furthermore, 5′ and 3′ RACEs (rapid amplification of cDNA ends) procedures were performed with two primers pairs (available upon requested) to obtain the full sequence as described previously (Lan and Lu, 2020). ORFs were predicted on the NCBI website, https://www.ncbi.nlm.nih.gov/orffinder/. Sequence identity comparison of HCRSV-ZZ 5′-terminal, seven ORFs and 3′terminal regions was performed with software BioEdit. The phylogenetic tree was established by maximum likelihood method with MEGA 5.1 software with 1,000 replicates. Prediction of hairpin loop 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 virus-infected leaves of H. rosa-sinensis plants were made on 200 mesh Formvar coated copper grids and stained with 5% (wt/vol) uranyl acetate. Virus particles were observed under H-7650 Hitachi transmission electron microscope (Hitachi, Tokyo, Japan) at 80 kV as described previously (Lan and Lu, 2020; Lan et al., 2019).

Total RNA was extracted from virus-infected and healthy leaves of hibiscus plants, H. 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 and Lu, 2020; Lan et al., 2015, 2016, 2019). Briefly, small RNA molecules (18-32 nt) were purified and Solexa adaptors were ligated to 5′- and 3′-terminals of small RNA molecules. Then, 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 (San Diego, CA, USA). 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 HCRSV genome. The downstream analyses of vsiRNAs were performed using Perl scripts and Excel.

H. rosa-sinensis plant samples showing severe chlorotic mottle and curl of plant leaves (Fig. 1A) for about half a month were collected from Taiwan Farmers Entrepreneurship Park in Zhangzhou city of Fujian province of China. HCRSV infection was confirmed by RT-PCR and H-7650 Hitachi transmission electron microscopy successively. It was firstly showed that CP gene (about 500 bp) of HCRSV was successfully amplified by RT-PCR from H. rosa-sinensis leaves showing severe chlorotic mottle and curl symptom but not from healthy plant leaves with specific primers (Fig. 1B). It was secondly suggested that about a large amount of virus particles with 30 nm in diameter were found in leaf dips prepared from H. rosa-sinensis leaves showing severe chlorotic mottle and curl symptom (Fig. 1C) but not from healthy plant leaves (data not shown) by H-7650 Hitachi transmission electron microscopy. Taken together, these results confirmed that H. rosa-sinensis plants were indeed infected by HCRSV pathogen.

The complete nucleotide sequence of HCRSV isolated from Zhangzhou city of Fujian province in China (named HCRSV-ZZ) and its deduced amino acid sequence were obtained. The complete nucleotide sequence of HCRSV-ZZ was 3,909 nt in length, including 5′-UTR, internal seven ORFs and 3′-UTR (Table 1, Fig. 2).

The length of the 5′-UTR of HCRSV-ZZ isolate was 30 nt, which is longer than that of HCRSV Colombia (19 nt), Brazil (23 nt), South Africa (19 nt), and Malaysia (18 nt) isolates and similar to that of HCRSV Xiamen (30 nt), Iran (30 nt), Italy (29 nt), Israel (30 nt), America (30 nt), Singapore (30 nt), and Taiwan (30 nt) isolates. The 5′-UTR of HCRSV-ZZ isolate shared the highest (94.4%) and the lowest (76.7%) identity to that of HCRSV Colombia and Taiwan isolates at the nucleotide level, respectively (Table 1). The 5′-UTR of HCRSV-ZZ isolate also contained well-characterized and conserved sequences (Supplementary Fig. 1), which has been shown to function as subgenomic promoters during regulating viral RNA synthesis (Li and Wong, 2006).

The genome of HCRSV-ZZ contained seven ORFs encoding for P28 (nt 31-675), P23 (nt 41-670), P81 (nt 31-2,238), P8 (nt 2,205-2,414), P9 (nt 2,341-2,577), P38 (nt 2,587-3,624), and P25(nt 2,600-3,274), respectively (Table 1, Fig. 2) by being predicted on the NCBI website. ORF(P23) was within ORF(P28) (Fig. 2). P81 was generated by in-frame readthrough of the stop codon of ORF(p28) towards the 3′-terminus of the viral genome (Fig. 2). ORF(P81) and ORF(P8) overlapped partially by 34 nt each other (Fig. 2). ORF(P8) and ORF(P9) overlapped partially by 74 nt each other and located at the central of the viral genome (Fig. 2). ORF(P38) was located at 3′-terminus of the viral genome (Fig. 2). ORF(P25) was within ORF(P38) (Fig. 2). All seven ORFs of HCRSV-ZZ isolate shared the highest identity to that of HCRSV Colombia isolate and the lowest identity to that of Taiwan isolates at the nucleotide and protein level, respectively (Table 1).

The length of the 3′-UTR of HCRSV-ZZ isolate was 285 nt (Table 1, Fig. 2), which is shorter than that of HCRSV Colombia (319 nt), Brazil (325 nt), and America isolates (342 nt) but longer than that of HCRSV South Africa (179 nt) and Malaysia (258 nt) and similar to that of HCRSV Xiamen (284 nt), Iran (284 nt), Italy (284 nt), Israel (284 nt), Singapore (284 nt), and Taiwan (283 nt) isolates. The 3′-UTR of HCRSV-ZZ isolate shared the highest (98.6%) and the lowest (85.4%) identity to that of HCRSV Colombia and Taiwan isolates at the nucleotide level, respectively (Table 1). It was shown that a well-characterized and highly conserved sequences, 5′-AATACCCTCACAACCGCCATAGTTGTGAGGCAGCCC-3′, was also found in the 3′-UTR of HCRSV-ZZ isolate as other members of the genus Betacarmovirus (Supplementary Fig. 2). The sequence could be folded into a small stable hairpin loop structure (Supplementary Fig. 2) by predicted with RNAfold server (http://rna.tbi.univie.ac.at/), which has been shown to play an important role in the replication of viral gRNA and sgRNAs (Li and Wong, 2006). Taken together, these results suggested that HCRSV-ZZ have a genome organization or arrangement similar to that of other viruses in genus Betacarmovirus.

Furthermore, it was shown that HCRSV-ZZ isolate clustered closely with HCRSV Colombia isolate (Fig. 3) by the maximum-likelihood phylogenetic trees inferred from the whole viral genome sequences, implying high homology among them.

To characterize the HCRSV-derived vsiRNAs from H. rosa-sinensis, small RNA libraries were constructed and sequenced with the Solexa protocol (Gen Denovo) with total RNA extracted from HCRSV-ZZ-infected and healthy H. rosa-sinensis. A total of 11,049,244 and 8,733,680 siRNAs reads were sequenced from HCRSV-ZZ-infected and healthy H. rosa-sinensis, respectively. It was shown that 24 nt siRNAs were most dominant in HCRSV-free H. rosa-sinensis, but 21 nt and 22 nt siRNAs were most dominant in HCRSV-infected H. rosa-sinensis (Fig. 4A), which is similar to previous reports for several pathogen systems (Lan and Lu, 2020; Lan et al., 2019). That is to say, the reads of 21 nt and 22 nt siRNAs increased, but 24 nt siRNA decreased in HCRSV-infected H. rosa-sinensis (Fig. 4A). These changes in most abundance of siRNAs showed that HCRSV infection could modulate the length distribution pattern of siRNAs in H. rosa-sinensis, or HCRSV infection triggered the antiviral immunity response in this plant.

By being mapped to the HCRSV genome, 11,040 unique vsiRNAs reads (18-32 nt, accounting for 0.10% of all unique siRNAs reads) were characterized from the HCRSV genome in H. rosa-sinensis infected by HCRSV (Fig. 4B). However, a very small number of 138 unique vsiRNAs reads (18-32 nt, accounting for 0.0016% of all unique siRNAs reads) were characterized from the HCRSV genome in healthy H. rosa-sinensis (data not shown). In HCRSV-infected H. rosa-sinensis, 21 nt, 22 nt, and 20 nt vsiRNAs were most abundance, accounting for about 48.3% of vsiRNAs reads (Fig. 4B). In numerous pathogen systems, the majority of vsiRNAs were 21 nt and 22 nt in length (Kreuze et al., 2009; Lan and Lu, 2020; 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, length distribution of 21 nt, 22 nt, and 20 nt vsiRNAs derived from HCRSV implied that the homologs of DCL4 and DCL2 proteins in H. rosa-sinensis may respectively be also the predominant DCLs involved in the biogenesis of these vsiRNAs which functions as the predominant antiviral silencing mediators as previous reports (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 H. rosa-sinensis participating in generation of vsiRNAs with different lengths will be the future discussed topic. Furthermore, the overhangs (No. of unpaired nucleotides at 5′ or 3′-terminal) of 21 nt vsiRNAs duplexes were analyzed. It was shown that 21 nt vsiRNA duplexes with 2 nt overhangs were the most abundant, followed by 21 nt duplexes with 1 nt overhangs in H. rosa-sinensis infected by HCRSV (Fig. 4C), which was similar that of other viruses-plant and viruses-invertebrates systems (Liu et al., 2018; Niu et al., 2017). Thus, these results revealed that 21 nt vsiRNAs duplexes with 2 nt overhangs, generated by the homologs of DCL4 and DCL2 proteases, were the most efficient triggers of RNAi by HCRSV infection in H. rosa-sinensis.

Previous reports have shown that the 5′-terminal nucleotides of vsiRNAs play an irreplaceable role in controlling the sorting of vsiRNAs to different AGO proteins in plants (Mi et al., 2008). In this study, it was shown that HCRSV vsiRNAs (18-32 nt) demonstrated a clear tendency to begin sequentially with uracil (U), cytosine (C), adenine (A), and guanidine (G) (Fig. 4D), which was consistent with previous studies for diverse plant-virus systems (Donaire et al., 2008, 2009; Lan and Lu, 2020; Lan et al., 2019; Xu and Zhou, 2017). To obtain further understanding of HCRSV vsiRNAs sorting, the 5′-terminal nucleotides of 21 nt, 22 nt, and 20 nt vsiRNAs, the top three abundant in reads, were furthermore analyzed. It was suggested that these vsiRNAs (21 nt, 22 nt, and 20 nt) showed a similar tendency to begin sequentially with U, C, A, and G (Fig. 4D). The preference the 5′-terminal nucleotides of vsiRNAs is consistent with the AGO1 function in defending against RNA viruses in plant (Morel et al., 2002; Qu et al., 2008). However, for 23 nt, 24 nt, and 25 nt HCRSV vsiRNAs, it was shown that a strong bias for vsiRNAs beginning with a 5′-A was observed (Fig. 4D), implying the high binding affinity of AGO2 and AGO4 with these vsiRNAs (Mi et al., 2008). Additionally, the low proportion of vsiRNAs beginning with G was also observed (Fig. 4D), consistent with previous reports of pathogen-systems (Donaire et al., 2009; Li et al., 2016; Mi et al., 2008; Mitter et al., 2013). Thus, the preference the 5′-terminal nucleotides of vsiRNAs suggested that diverse AGO proteins were involved in vsiRNAs binding in H. rosa-sinensis infected by HCRSV. In the future, isolation and determination of AGO-vsiRNAs complexes from H. rosa-sinensis infected with HCRSV will provide us more detail functions information about vsiRNAs sorting in hibiscus plants.

To reveal the origin of vsiRNAs, vsiRNAs were located to positive and negative strands of the HCRSV genome. It was shown that 72% and 28% of HCRSV vsiRNAs were derived from the viral positive and negative strands of the HCRSV genome, respectively (Fig. 5A). The origin distribution pattern was consistent with that of hibiscus latent fort pierce virus (HLFPV)-derived vsiRNAs in the same plant host, H. rosa-sinensis, in which, HLFPV vsiRNAs also derived predominantly from the viral positive-strand RNA (Lan and Lu, 2020), implying the two viruses were encountered with the similar RNAi antiviral immunity mechanism in the same plant host, H. rosa-sinensis. Additionally, single-base resolution maps of total unique vsiRNAs along with the HCRSV genomes were established with Bowtie tools and in-house Perl scripts. It was shown that HCRSV vsiRNAs were derived from the whole HCRSV genome (Fig. 5B) and presented a discontinuous and heterogeneous distribution pattern (Fig. 5B), i.e., forming the high-frequency cutting site (Hot spot) and low-frequency cutting site (Cold spot) along the HCRSV genome (Fig. 5B). Taken together, the origin distribution pattern and single-base resolution maps implied that the regions of stem-loop secondary structures present in HCRSV genome may be the substrates for DCL proteins for generating vsiRNAs. To confirm the speculation, we evaluated the stem-loop secondary structures of HCRSV genome with RNAfold server. It was shown that lots of stem-loop secondary structures were formed in the whole genome regions and these structures corresponded to the hotspots (Fig. 5B-D). Taken together, these results suggested that HCRSV vsiRNAs should originate predominantly by direct DCL cleavage of imperfect duplexes in stem-loop secondary structures of the positive strand of HCRSV 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).