Highly divergent isolates of chrysanthemum virus B and chrysanthemum virus R infecting chrysanthemum in Russia

Sampling

We gathered leaf samples in the chrysanthemum (Chrysanthemum morifolium Ramat.) germplasm collection of the Nikita Botanical Gardens (44.51N; 34.23E) from August to October 2019–2020. Plants displaying symptoms of supposed virus infections were initially selected. Samples of plants with unclear symptoms and symptomless plants were also taken. Three to four leaves collected from each plant were pooled for subsequent analyses. Each plant was analyzed separately.

RNA isolation, RT-PCR, and Sanger sequencing

We tested leaf samples for CVR and CVB using RT-PCR. Total RNA was isolated using RNeasy Plant Mini Kit (Qiagen, Hilden, Germany). Random hexamer primers and MMLV reverse transcriptase (Evrogen, Moscow, Russia) were employed for the first-strand cDNA synthesis. To detect CVB, primers P1/P2 amplifying the entire CP gene were used for PCR according to the original protocol (Ram et al., 2005). CVR was detected using primers CV1zF/CV807R, which target the 5′-terminal sequence of the replicase gene (Wang et al., 2018). The primer CV1zF (5′-ATGGCGCTCACTTTCCGCAGCC-3′) is a modified version of the original forward primer CV1F (Wang et al., 2018) that was found to be fitted better for the successful detection of Russian CVR isolates. The HVR of the CVR replicase gene was amplified using the primers CV1406F/CV2353R as described previously (Wang et al., 2018). The 3′-terminal CVR genome region covering the entire TGB3 and CP genes was amplified using the forward primer CV7162F (Wang et al., 2018) and a newly designed reverse primer CVR-R (5′-TATGACCCCTAGCCTTTTGA-3′), which target ORF3 and ORF6 of the CVR genome, respectively. A number of primers, recognized different regions of the CVR replicase gene, were also designed to confirm the HTS results. Information about primers used in this work is summarized in Table S1. Amplification was performed using the proof-reading Encyclo DNA polymerase (Evrogen, Moscow, Russia) and included the following cycles: denaturation at 94 °C for 30 s, primer annealing at 56 °C for 30 s, and elongation at 72 °C for 1 to 2 min depending on the size of expected PCR product, for 35 cycles with a final extension at 72 °C for 10 min. Amplicons were analyzed by 1.5% agarose gel electrophoresis using ethidium bromide staining and MultiDoc-It Digital Imaging System (Analytik Jena GmbH, Jena, Germany). Appropriate PCR products were purified from agarose gel with BC022 Cleanup Standard Kit (Evrogen, Moscow, Russia) and were subjected to Sanger sequencing in both directions using Evrogen facilities. If several plants of the same cultivar were found to be infected, PCR products from only one of them were sequenced.

Assembly of viral genome sequences from RNA-Seq data

We generated three DNA libraries using polyA-enriched fraction of total RNA extracted from plants of the cultivars Ribonette, Fiji Yellow, and Golden Standard with mixed CVB and CVR infection. This fraction was obtained using a K0600 full-size poly(A) mRNA purification magnetic beads kit (Sileks, Moscow, Russia) according to the manufacturer’s protocol. TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA) was employed for the сDNA libraries synthesis, according to the manufacturer’s instructions. The libraries were sequenced on MiSeq Illumina platform. Reads were trimmed and filtered using Trimmomatic v.0.39 (Bolger, Lohse & Usadel, 2014). Next, contigs were assembled de novo using Trinity 2.11 with default settings (Grabherr et al., 2011). The virus-specific contigs were identified by BLASTn search using GenBank full-length viral genome database (https://ftp.ncbi.nlm.nih.gov/genomes/genbank/viral/assembly_summary.txt).

Analysis of assembled genomes

Multiple alignments of genomic sequences using ClustalW, determination of their nucleotide (nt) and amino acid (aa) identities as well as phylogenetic analysis using the neighbor-joining algorithm and Kimura 2-parametric model were performed by the MEGA7 program (Kumar, Stecher & Tamura, 2016). Sequences of all available CVR and CVB isolates were retrieved from GenBank and employed for the analyses. Recombination detection program RDP4.46 (Martin et al., 2015) was used to search for recombination in CVR and CVB genomes using default settings with the exception that ‘sequences are linear’ and ‘list events detected by > 6 methods’ options were chosen.

Prevalence of CVB and CVR in chrysanthemum collection

A total of 90 chrysanthemum plants of 22 most popular cultivars and one hybrid were surveyed (Table 1). By RT-PCR, CVB was detected in plants of 16 cultivars and hybrid 7–15. CVR was identified in eight cultivars and hybrid 7–15. PCR products of the expected sizes 948 base pairs (bp) and 747-bp were amplified in the samples infected with CVB and CVR, respectively (Fig. S1). A total of 52 plants (58%) were infected with CVB and 39 plants (43%) with CVR. The data presented in Table 1 show that CVB and CVR are mostly distributed across the collection independently to one another, although mixed infection was detected in 21 plants (23%). These viruses were detected in plants exhibiting mild mosaic, mottle or pale green spots and in symptomless plants as well suggesting that such the symptoms were non-specific for CVR and CVB (Fig. S2). In addition, these viruses were not detected in plants of the cultivar Yunost with symptoms of mottle on the leaves. No stunting was observed among tested plants. It is possible that the symptoms manifested in a cultivar-dependent manner. Neither CVB nor CVR were detected in the cultivars Dayon, Milka Pink, Plamya, and Yunost. Thus, we detected CVR in Russia and outside China for the first time expanding the information on the geographic distribution of the virus. We also demonstrated that CVR can be widespread on chrysanthemum.

Analysis of partial CVR genome sequences

Sequencing of 747-bp PCR products, generated in CVR-positive samples, showed that CVR isolates from the cultivars Sirenevye Dali, Skazka, Zolotovoloska, Golden Standard, Chita, Yuzhnaya Notch, and hybrid 7–15 were, in average, 96% identical to each other and the corresponding genomic region of the Chinese isolate BJ (MG432107). Until recently, the latter one was the only CVR isolate, which a full-length genome was available. At the same time, isolates from the cultivars Ribonette and Fiji Yellow shared 95% identity to each other and were only 83% identical to the rest ones and the isolate BJ, suggesting that they may represent a genetically distant CVR group. To test this assumption, we attempted to sequence two other genome regions of the Russian CVR isolates. The PCR products of the expected sizes 946 and 1,420-bp were generated by amplification of the hypervariable and 3′-terminal genomic regions, respectively, in most samples except the cultivars Fiji Yellow and Ribonette (Fig. S3). Partial genome sequences of CVR isolates from the cultivars Sirenevye Dali, Skazka, Zolotovoloska, Chita, Yuzhnaya Notch, and hybrid 7–15 were deposited in GenBank under accession numbers MZ231086 to MZ231090 and MZ340424 to MZ340436. We used these sequences for the phylogenetic analysis of three genomic regions. The HVR and TGB3-IGR2-CP sequences of the isolates from the cultivars Ribonette and Fiji Yellow were retrieved from their full-length genomes (see below). The trees were congruent and the isolates from the Ribonette and Fiji Yellow cultivars always formed a separate clade with 100% bootstrap values further supporting the idea that they belong to a divergent group of CVR (Fig. S4).

Analysis of partial CVB genome sequences

Sequences of 948-bp PCR products, generated in CVB-positive samples and encompassed the entire CP gene, were deposited in GenBank under accession numbers MZ231081–MZ231085. Their phylogenetic analysis together with the CP of isolates from Russia (MH678697−MH697704) (Zakubanskiy et al., 2021), India, Korea, China and Japan available in GenBank showed that most Russian isolates were distributed among three groups (A1, A2, and B) as supported by high bootstrap values (Fig. 1). Within a group, the nt identities of the Russian isolates ranged from 95.6% to 99.7% in group A1 and 93.5% to 98.5% in group A2. The CVB isolates from the cultivars Perfection Red and Gracia (group B) were 82.5% identical to one another and 77.5% to 83.6% to the isolates from groups A1 and A2, indicating higher CP variability in group B. The isolate from the cultivar Showdown White was out of these groups and was most closely related to a Chinese isolate (JQ904593), sharing 96% identity. We failed to sequence the CP gene of the isolate from the cultivar Golden Standard due to a large number of ambiguous positions at the corresponding chromatograms suggesting that this cultivar can be co-infected with divergent CVB isolates (Data S1).

Metatranscriptomic sequencing of samples from the cultivars Ribonette, Fiji Yellow, and Golden Standard

Partial genome sequences analysis showed that CVR isolates from the cultivars Ribonette and Fiji Yellow (named X6 and X13, respectively) appear to represent a divergent group of the virus, while that from the cultivar Golden Standard (isolate X21) is closer to other Russian isolates and the Chinese isolate BJ. On the other hand, the CVB isolates from the cultivars Ribonette and Fiji Yellow (named Rib and FY, respectively) are assigned to different phylogenetic groups considering the CP sequences. The cultivar Golden Standard is probably infected with divergent CVB isolates. Based on partial genome sequences, the samples from the cultivars Fiji Yellow, Golden Standard, and Ribonette were selected to determine CVR and CVB complete genomes using HTS.

A total of 12.4, 13.7, and 25.6 million of pair-end 150 nt reads were generated in the samples from the cultivars Ribonette, Fiji Yellow, and Golden Standard, respectively. CVB-, CVR-, and TAV-related reads were identified in each of these three samples. The proportion of virus-related reads is presented in Table 2. The raw reads were deposited in the NCBI Sequence Read Archive (SRA) (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA751454/). No reads related to other viruses were detected by RNA-Seq analysis. Further we focused on the characterization of two carlaviruses genomes.

The assembled CVR-related contigs covered the whole genome of the reference isolate BJ (MG432107, 8,990 nt) except 4 to 12 nt at the 5′-end of the genome. Five genome regions of the total length 4,046 nt (about of half of the genome) were re-sequenced by the Sanger method using primers presented in Table S1. The sequences determined by HTS and the Sanger method were shown to be identical. The near complete genomes of CVR isolates from the cultivars Ribonette (isolate X6; 8,878 nt), Fiji Yellow (X13; 8,872 nt), and Golden Standard (X21; 8,873 nt) were deposited in GenBank under accession numbers MZ514905 to MZ514907.

The assembled CVB-related contigs from the cultivars Golden Standard and Fiji Yellow covered the whole genome of the reference isolate S (AB245142, 8,990 nt) except 12 to 19 nt at the 5′-end of the genome. Two different CVB genomes (designated GS1 and GS2) were assembled in the cultivar Golden Standard. The near complete genomes of the CVB isolates GS1 (8,977 nt), GS2 (8,976 nt), and FY (8,970 nt) were deposited in GenBank under accession numbers MZ514908 to MZ514910. The CVB genome from the cultivar Ribonette (isolate Rib) was failed to assemble completely. The partial genomic sequence of 5,174 nt covering the C-terminal part of the replicase gene, TGB, CP, and CRP genes as well as 3′-UTR was deposited in GenBank under accession number MZ514911. The CP sequences of the isolates FY and Rib determined by the HTS and the Sanger method were identical.

Characterization of new CVR genomes

The structure of newly sequenced CVR genomes was typical of carlaviruses. These comprise six ORFs, 5′- and 3′-UTRs, and two IGRs between ORF1 and ORF2 as well as between ORF4 and ORF5. Complete genomes of the Russian isolates and Chinese isolate BJ were compared (Table 3). Isolate X21 was 92% identical to the BJ and similar to the latter in the size and organization of all genomic elements. In contrast, genomes of the isolates X6 and X13 were 76% to 77% identical to the X21 and BJ and shared 95% identity to one another. ORF4 was the most variable genome region both at the nt and aa levels. In opposite, ORF5 and ORF6 were the most conserved. Phylogenetic analysis showed that isolates BJ, X21 and X6, X13 formed two separate clades with 100% bootstrap support (Fig. 2). Using NCBI Conserved Domain Search Service (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi), the MTR, OXY, P-Pro, HEL and RdRp domains were identified in the X6, X13, and X21 replicase proteins at identical or close positions. At the same time, multiple alignment of ORF1s revealed a number of indels that resulted in several aa deletions located between MTR, OXY, and P-Pro domains. Another aa deletion was predicted in HEL domain of the X6 and X13 (Fig. S5). These indels were confirmed by Sanger sequencing of the corresponding genome regions using appropriate primers (Table S1). Due to these distinctions, X6 and X13 ORF1s are shorter by 6 nt and the encoded replicase protein is shorter by 2 aa when comparing with the isolates X21 and BJ. Furthermore, IGR1s of the X6 and X13 are of 27 nt length that shorter by 1 nt than those in the isolates X21 and BJ. Also, ORF3 and ORF5 in the isolates X21 and BJ are terminated with TAA and OFR6 with TGA. On the contrary, in the X6 and X13 the stop codon in ORF3 and ORF5 is TGA, whereas ORF6 is terminated with TAA. Because ORF5 and ORF6 overlap, replacing the TAA terminating codon with TGA in ORF5 of the X6 and X13 shifts putative ORF6 start AUG codon by 3 nt upstream. This shift resulted in the CRP elongation by 1 aa comparing to the X21 and BJ and in changes of aa sequence in the N-terminal region of the protein (Fig. 3). However, neither NLS (⁴⁸RRRR⁵¹) nor ZF motifs (⁵⁸CX₂CX₁₂CX₄C⁷⁹) were affected in the isolates X6 and X13.

Thus, based on near complete genome sequences and phylogenetic analysis, at least two divergent groups of CVR exist, which can be probably considered different strains of the virus.

In addition, the X6 and X13 genome analysis showed that these isolates were obviously not recognized by primers targeting HVR and 3′-terminal genome regions due to a large number of mismatches between the primers and the corresponding genome sequences (3 to 9 per a primer). In contrast, RT-PCR with primers CV1zF/CV807R, specific to the beginning of the polymerase gene, recognized even highly divergent isolates of the virus and apparently can be successfully applied for the CVR detection.

Characterization of new CVB genomes

We assembled near complete genome sequences of the CVB isolates FY, GS1, and GS2. Typically for carlaviruses, six ORFs, 5′- and 3′-UTRs, and three IGRs between ORF1 and ORF2, ORF4 and ORF5, ORF5 and ORF6 were identified in the new genomes.

Seven complete CVB genomes were available in GenBank until recently. Four of these are from India (AM493895, AM765837−AM765839) (Singh et al., 2012), two were found in China on a new host Gynura japonica (MW269552–MW269553) (Lai et al., 2021), and another isolate S is from Japan (AB245142) (Ohkawa et al., 2007). The multiple alignments showed that the isolates GS1 and GS2 were 82% identical to each other and shared 82% to 84% identity with the FY. Overall, the Russian CVBs were 83% to 87% identical to the Chinese and Japanese isolates and 68% to 72% identical to the Indian ones. The Japanese isolate S was used as a reference to compare the Russian complete genomes (Table 4). ORF1 is the most variably part of the FY and GS2 genomes while ORF6 of the FY and GS1 is the most conserved. At the same time, ORF4 of the FY and ORF6 of the GS2 were only distantly related to their counterparts from other isolates. Phylogenetic analysis of the complete or near complete genomes of all available CVB isolates showed that the Indian isolates form a separate cluster (Fig. 2). Another cluster is represented by isolates from Russia, Japan, and China. The FY and GS2 form a separate clade, and the isolate GS1 is clustered with the Japanese isolate S.

In addition, the detection of two different CVB isolates, GS1 and GS2, in the same plant of the cultivar Golden Standard by HTS enables to explain why we failed to obtain a consensus sequence of the CP gene by Sanger sequencing. A large number of ambiguous positions at chromatograms is probably because the CPs of the GS1 and GS2 are identical only 83% at the nt level.

Several recombination events (RE) were detected in the genomes of Indian and Japanese isolates previously (Singh et al., 2012). In this work, we looked for recombination in the extended alignment that included near complete genomes of new Chinese and Russian CVB. Three REs were predicted in the replicase gene (Fig. 4). The RE1 in the FY genome was supported by seven independent algorithms implemented in RDP4 with statistically significant P-values (RDP = 3.4 × 10⁻²⁰⁸, GENECONV = 4.9 × 10⁻¹⁶⁹, BOOTSCAN = 7.2 × 10⁻¹⁹⁰, MAXIMUM CHI SQUARE = 1.5 × 10⁻⁴⁹, CHIMAERA = 5.4 × 10⁻⁴⁹, SISCAN = 7.1 × 10⁻¹²¹, 3SEQ = 10⁻¹⁴⁶, and 2.9 × 10⁻¹⁶⁹ in average). The RE1 was bounded with breakpoints at positions 602 and 4,322, encompassing C-terminal part of MTR domain, the entire OXY and P-Pro domains, and the N-terminal part of HEL domain. The recombinant sequence is 3,720 nt, which is about 40% of the whole genome. The inferred major parent was the Indian isolate Uttar Pradesh (AM765837), and the minor parent was unknown. However, the isolates FY and Uttar Pradesh have only 88% identity beyond the recombinant part of their genomes. In addition, a BLASTn search along GenBank nucleotide collection (https://blast.ncbi.nlm.nih.gov/Blast.cgi) showed that the recombinant sequence was most closely related to the corresponding genomic region of the CVB isolates HZ-V2 sharing 79% identity. These data suggest that some other CVB isolates seemed to be the parents of the recombinant isolate FY. Apparently, closer isolates have to be found yet or their full-length genomes are not available at the moment. The RE1 was also confirmed by phylogenetic analysis of CVB isolates (Fig. S6). The trees reconstructed from different genome regions corresponding to the putative recombinant sequence of the isolate FY and to the rest of the genome were incongruent suggesting recombination (Chare & Holmes, 2006). The RE2 and RE3 were also identified in the isolates Punjab and HZ-V2 with a high degree of confidence (the average P-values were 3.8 × 10⁻¹⁵¹ and 1.2 × 10⁻¹³⁰, respectively).