Genetic and phylogenetic analysis of Chinese sacbrood virus isolates from Apis mellifera

Clinical sample collection and testing

In total, 359 samples (each collected sample included five larvae from a single colony) were obtained from China between 2008 and 2017. Among all 359 samples, 20 samples were from A. mellifera and 359 samples were from A. cerana. In all larvae, CSBV infection was detected by real-time quantitative RT-PCR (qRT-PCR) (Ma et al., 2013a). The field studies did not involve endangered or protected species, and the owner of Linyi bee farm Zhang Youcai gave permission to conduct the study on this site.

We randomly selected one positive larva for genetic characterization of VP1 from the same colony at the same time. The VP1 gene of CSBV was amplified using VP1-specific primers (F: 5′-GCGGATCCATGGATAAACCGAAGGATATAAG-3′, R: 5′-GCAAGCTTTTATTGTACGCGCGGTAAATA-3′) and sequenced.

Phylogenetic tree construction from VP1

Using VP1 as a target gene, we constructed a phylogenetic tree from all the CSBV isolates and reference strains in GenBank (Tables 1 and 2) using the MEGA 5.0 package (Tamura et al., 2007) and the neighbor-joining (NJ) method (Saitou & Nei, 1987). The phylogenetic tree was computed using the Kimura 2-parameter method (Kimura, 1980). A bootstrap value of 1,000 replicates was applied to yield a robust phylogeny. We named the samples after the strain that originated in the same region and showed 100% homology, then submitted the data to GenBank (CSBV uniformly renamed as CSBV + first isolated geographic location + year).

Virus purification

AmCSBV-SDLY-2016 was obtained from a natural outbreak on Linyi bee farm in Shandong, China. A. mellifera larvae infected by AmCSBV-SDLY-2016 were collected by the owner of Linyi bee farm. A total of 50 infected A. mellifera larvae displaying upward warping of the head and body surface changes, and white or sick and dead larvae in the honeycomb with the capped brood abnormally sealed, were collected, weighed, and completely homogenized in sterile water (1.5× by weight) using a mortar and pestle. AmCSBV-SDLY-2016 purification was performed by cesium chloride gradient centrifugation, according to the method reported by Ma et al. (2011b, 2011c). The supernatant was extracted with an equal volume of 1,1,2-trichlorotrifluoroethane before the aqueous phase was layered over a discontinuous CsCl gradient (1.5 and 1.2 g/cm³) and centrifuged at 270,000×g for 1 h. The material at the CsCl interface was harvested and adjusted to a volume of five mL (final density 1.38 g/cm³) with CsCl solution and centrifuged at 270,000×g overnight. The supernatant was then successively passed through 0.45 and 0.22-μm cell filters. Healthy larvae treated by the same method were used as the negative control. CSBV was subsequently identified by reverse-transcription polymerase chain reaction (RT-PCR) to exclude black queen cell virus (BQCV), acute bee paralysis virus (ABPV), chronic bee paralysis virus (CBPV), deformed wing virus (DWV), Kashmir bee virus (KBV), and Israeli acute paralysis virus (IAPV), following the method reported by Yu, Liu & Wang (2016). CSBV virus samples free of other viruses were stored at −80 °C until further use.

Electron microscopy and sodium dodecyl sulfate polyacrylamide gel electrophoresis for virus identification

As previously described (Ma et al., 2011b, 2011c), 100 μL of the purified viral suspension was directly pelleted onto carbon-coated Formvar copper grids by ultracentrifugation (15 min at 82,000×g) using a Beckman Airfuge. The grids were negatively stained with 2% sodium phosphotungstate at pH 6.8 for 90 s and observed using a Philips CM10 transmission electron microscope.

Structural proteins were separated by sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) with 5% stacking and 12% separating gels using standard protocols.

Agar gel immunodiffusion assay

Briefly, one g of agarose and eight g of sodium chloride (NaCl) were added to 100 mL phosphate buffer (0.01M, pH 7.2), shaken well, and microwaved for 2 min to prepare an agar solution. The solution was slightly cooled, poured into Petri dishes (90 mm in diameter; 20–22 mL of agar per plate), and allowed to solidify. Seven wells were made in the agar plates to identify AmCSBV-SDLY-2016. The central hole was loaded with antisera against the CSBV-JLCBS-2014 strain; the surrounding wells 1 and 2 were loaded with purified AmCSBV-SDLY-2016; wells 3 and 4 were loaded with a known positive control (purified CSBV-JLCBS-2014); and wells 5 and 6 were loaded with treated healthy larvae as the negative control. The plates were then placed in a closed, moist container and incubated at 37 °C in a humidified chamber for 24 or 48 h. Precipitation was visible after 24 h but became more distinct after 48 h of incubation.

Analysis of pathogenicity

As described previously (Ma et al., 2013a), PCR probes were used to measure the copy numbers of AmCSBV-SDLY-2016. The purified AmCSBV-SDLY-2016 was 10-fold serially diluted.

As described previously (Hu et al., 2016), two-day-old A. cerana larvae from a single mated queen honeybee in a healthy apiary in Jinzhou, Liaoning Province, were orally inoculated with CSBV. In total, 120 larvae were selected and randomly distributed into six groups (n = 20 larvae per group). Each group was inoculated with AmCSBV-SDLY-2016 (1.25 × 10⁴, 1.25 × 10⁵, 1.25 × 10⁶, 1.25 × 10⁷, and 1.25 × 10⁸ copies/larva in groups 1–6, respectively). The assay was performed in triplicate. The virus inoculation assay was performed twice after the end of the first experiment.

Each larva was fed 20 μL of a virus suspension mixed with an equal amount of basic larval diet (BLD, Liu & Zeng, 2010) consisting of 50% royal jelly, 37% sterile water, 6% glucose, 6% fructose, and 1% yeast extract, at 95% relative humidity and 34 °C. The virus-free control was fed 20 μL of sterile water with an equal amount of BLD. BLD was used subsequently for daily feeding. The clinical signs in each group of larvae were examined and recorded every day until larvae death. All larvae that were orally inoculated in this study were analyzed by RT-PCR for BQCV (Grabensteiner et al., 2007), ABPV (Grabensteiner et al., 2007), CBPV (Tentcheva et al., 2004), DWV (Tentcheva et al., 2004), KBV (Blanchard et al., 2014), IAPV (De Miranda, Cordoni & Budge, 2010), and qRT-PCR for CSBV (Ma et al., 2013a).

The statistical analysis of larval mortality was performed between groups and between three replicates of each group using SPSS 22.0.

AmCSBV-SDLY-2016 genome sequencing

The primers used in this study (Table 3) were designed based on the nucleotide sequences of CSBV-JLCBS-2014, AmSBV-UK-2000, and CSBV-GD-2002. The full-length AmCSBV-SDLY-2016 genome was determined by Sanger sequencing of the RT-PCR fragments and 3′ rapid amplification of cDNA ends (RACE) (Clontech-Takara, Mountain View, CA, USA), according to the method described by Ma et al. (2011c). The PCR amplification product was cloned into the pMD-18-T vector (Takara Biotechnology Co. Ltd., Dalian, China). The plasmids were then used to transform Escherichia coli DH5α cells (Takara Biotechnology Co. Ltd., Dalian, China). For each RT-PCR fragment, five clones were randomly chosen and sequenced by Sangon Biotech Co. Ltd. (Shanghai, China). If all five clones were identical, the sequence was considered correct. The correct nucleotide sequences from all the fragments were assembled to build a continuous complete sequence using the DNASTAR software, and the VP1 gene sequence was analyzed as described above.

Phylogenetic tree construction from AmCSBV-SDLY-2016 and SBV/CSBV genome sequences

Nucleotide sequences of the amplified RT-PCR fragments were assembled to generate the entire genome of AmCSBV-SDLY-2016 using the DNASTAR program. Multiple nucleotide and deduced amino acid sequence alignments identified sequence variations between clones covering the same SBV genome regions using ClustalW in the MegAlign program (DNAStar Inc., Madison, WI, USA) and the published SBV/CSBV sequences. A phylogenetic tree was constructed from the nucleotide sequences of the coding regions of 43 previously reported SBV/CSBV strains from various countries and the AmCSBV-SDLY-2016 isolated in this study. The phylogenetic tree constructed using the MEGA 6.0 package (Tamura et al., 2007) and NJ method (Saitou & Nei, 1987) was computed using the Kimura 2 parameter method (Kimura, 1980). The phylogenetic tree was bootstrapped 1,000 times and bootstrap values placed on each branch.

Similarity and virus recombination analysis

Based on the results of the phylogenetic analysis, the isolated CSBV/SBV strains were divided into nine groups as follows: A (CSBV-FZ-2014), B (AmCSBV-SDLY-2016), C (CSBV-SXnor-2012), D (CSBV-SXYL-2015), E (CSBV-BJ-2012), F (CSBV-LNQY-2008), G (CSBV-JLCBS-2014), H (CSBV-GD-2002), and I (AmSBV-UK-2000) for similarity analysis and virus recombination analysis. The Recombination Detection Program (RDP) values were computed using the RDP 4.0 software for the preliminary screening of the major parent and minor parent sequences. Similarity plots and BootScan were computed using Simplot software (Lole et al., 1999) with the following parameters: a window of 200 base pairs (bp; step: 20 bp), with gap-stripping and Kimura (2-parameter) correction, using AmCSBV-SDLY-2016, CSBV-JLCBS-2014, CSBV-GD-2002, and AmSBV-UK-2000 for Find Sites, and AmCSBV-SDLY-2016 and CSBV-JLCBS-2014 as the query sequences.

Characteristics of the VP1 gene

Of the 359 samples tested, 326 were found to be positive by strand-specific qPCR, 18 of which originated from A. cerana, whose infection rate was 90.0% (18/20), and the 95% confidence intervals for infection rate were 69–99%, with 308 positive samples originating from A. mellifera whose infection rate was 90.9% (308/339), and the 95% confidence intervals for the infection rate were 85.90–95.81%. Positive samples originating in the same region that showed 100% homology were defined as a strain. A total of 31 CSBV strains were isolated according to the aforementioned virus-naming scheme.

The phylogenetic tree of VP1 revealed two clusters (Fig. 1). One was related to SBV strains that originated from A. cerana (AC genotype), the other was related to the SBV strains that originated from A. mellifera (AM genotype). The phylogenetic tree also showed that AmCSBV-SDLY-2016 originating from A. mellifera belonged to the clade containing the CSBV and other Asian strains. To simplify the naming of the strain, we used AmSBV and AmCSBV to represent the SBVs/CSBV strains isolated from A. mellifera; all the others represent the SBVs/CSBV strains that were isolated from A. cerana. The corresponding host and gene numbers are listed in Tables 1 and 2.

AmCSBV-SDLY-2016 identified by electron microscopy and SDS-PAGE

Electron microscopy showed large amounts of typical CSBV particles in the preparations from virus-infected larvae; CSBV particles were icosahedrons and had an approximate diameter of 26 nm (Fig. 2). No virus particles were observed in the control preparations from healthy larvae. The four main structural proteins of CSBV were separated by SDS-PAGE (Fig. 3). The molecular weights of the four proteins were about 44.2, 37.8, 31.5, and 30.5 kDa, respectively (Feng et al., 1998; Ma et al., 2011a).

AmCSBV-SDLY-2016 identified by AGID assay

Agar gel immunodiffusion (AGID) assay revealed the distinct CSBV-specific lines of precipitin observed between the wells containing AmCSBV-SDLY-2016 and the antisera against CSBV-JLCBS-2014, as well as the positive control wells (Fig. 4). No precipitin lines were observed in the negative control wells.

Pathogenicity

After oral inoculation with AmCSBV-SDLY-2016 (groups 1–5), all larvae were analyzed by RT-PCR and qRT-PCR. All of the other honeybee viruses were undetectable, whereas CSBVs were detectable. All the larvae in the virus-free control (group 6) lacked common honeybee viruses.

Three repeated experiments showed that upon inoculation with 1.25 × 10⁴, 1.25 × 10⁵, or 1.25 × 10⁶ copies of AmCSBV-SDLY-2016, larvae mortality rates were 35–45%, 65–75%, and 80–90%, whose 95% confidence intervals for the proportion dead pupae were 19–64%, 46–88%, and 62–97%, respectively (Table 4). In contrast, mortality of the larvae was 100% in the two groups, where each larva was sequentially inoculated with 1.25 × 10⁷ and 1.25 × 10⁸ copies of AmCSBV-SDLY-2016, whose 95% confidence intervals for the proportion dead pupae were 83–100%. In the virus-free control groups, the mortality rates of the larvae were 20–25%, with 95% confidence intervals for the proportion dead pupae of 6–44%.

Statistical analysis showed that there were no significant differences between the larval mortality from three repeated experiments in each group, but there were highly significant differences in larval mortality between each group, except that there was no difference in the number of larval mortality between the two groups, where each larva was sequentially inoculated with 1.25 × 10⁷ and 1.25 × 10⁸ copies of AmCSBV-SDLY-2016.

Nucleotide sequence

The nucleotide sequence of the AmCSBV-SDLY-2016 genome was 8,794-bp long, and the percentages of A, U, G, and C were 29.95%, 29.22%, 24.30%, and 16.52%, respectively, similar to those reported for other CSBV/SBV strains encoding one large ORF. Two AUG codons were located at positions 189 and 429 of the AmCSBV-SDLY-2016 genome; however, AUG 189 is likely the translation initiation site, as unlike AUG 429, it was observed in a sequence (AUUAUGG) identical to that of many invertebrate initiation codons (ANNAUGG). The size of the CSBV 3′ untranslated region (UTR) (77 nucleotides) was similar to that from other picornaviruses (40–126 nucleotides). The 5′ UTR was generated from AmCSBV-SDLY-2016 by RT-PCR with a primer designed using the CSBV-JLCBS-2014 sequence, indicating that the CSBV 5′ UTR was similar to that of CSBV-JLCBS-2014 (188 nucleotides). Multiple sequence comparisons revealed a sequence homology of 92.4– 97.1% among all CSBV isolates and a similarity of 94.5–97.7% in the deduced amino acid sequences. AmCSBV-SDLY-2016 was least similar (89.5–90.4% identity) to other SBVs but showed maximum similarity with CSBV-FZ-2014 (97.1% and 97.7% homology of sequence and deduced amino acid sequences, respectively). The results of VP1 sequencing were completely consistent, so we could rule out the possibility that several CSBV isolates could co-exist in the samples originating from A. mellifera.

Protein sequence

The deduced amino acid sequences of AmCSBV-SDLY-2016 genomes and previously reported SBV/CSBV strains were aligned and compared. Results revealed that the structural and non-structural proteins were located at the 5′ and 3′ ends, respectively. Multiple sequence alignment showed that AmCSBV-SDLY-2016 had 17 and one amino acid deletions (positions 711–729, which are in the region of structural protein VP1, and 2,128, which is in the region of nonstructural protein) as compared with CSBV-GD-2002, and three and 13 amino acid deletions (positions 711–713 and 715–728, which are in the region of structural protein, respectively) as compared to AmSBV-UK-2000 (Fig. 5). However, AmCSBV-SDLY-2016 was similar to the CSBV-JLCBS-2014 strain that infects A. cerana.

The amino acid sequence at the C-terminal region of the AmCSBV-SDLY-2016 polyprotein was similar to that for helicase, protease, and RNA-dependent RNA polymerase (RdRp) domains of the previously reported SBV/CSBV strains. The highly conserved amino acid motifs GPAGIGKS, QPVVVYDD, and KKIRGNPLIVILLCNH, corresponding to the helicase domains A, B, and C (Supplemental Information 1), respectively, were located between the amino acid positions 1,353 and 1,473; however, the C domain containing KKIRGNPLIVILLCNH appeared to be the most conserved, unlike mammalian picornaviruses. The conserved cysteine protease motif ²²⁴⁹GXCG²²⁵² (GACG) and the putative substrate-binding residue ²²⁶⁶GxHxxG²²⁷¹ (GMHFAG) were detected in the 3C protease domains spanning amino acids 2,141–2,272 (Supplemental Information 2). In addition, eight conserved domains identified in RdRp were also detected between amino acid positions 2,444 and 2,813 in AmCSBV-SDLY-2016 (Supplemental Information 3).

Phylogenetics

To assess the genetic relationship among the SBV/CSBV strains, a phylogenetic tree based on the nucleotide sequence of the SBV/CSBV-coding region was constructed by the NJ method for AmCSBV-SDLY-2016 and all SBV/CSBV reference strains in GenBank. The tree was similar to the phylogenetic tree of VP1, which revealed two clusters, one related to SBV strains that originated from A. cerana (AC genotype), and the other related to the SBV strains that originated from A. mellifera (AM genotype) (Fig. 6). The AM cluster was distinctly separate from the AC lineages, as supported by a clade credibility value of 100.0. The AC genotype was further subdivided into several subtypes according to their countries of origin and host species.

The phylogenetic tree also showed that AmCSBV-SDLY-2016 belonged to the clade containing the CSBV and other Asian strains. These results indicate that the AmCSBV-SDLY-2016 strain showed close genetic relationships with CSBV strains, specifically CSBV-FZ-2014.

Similarity profile

To identify differences in the full-length sequences and various genomic regions from CSBV strains AmSBV-UK-2000 and AmCSBV-SDLYA-2016, the complete coding regions were plotted using SimPlot, with AmCSBV-SDLY-2016 as the query sequence (Fig. 7). AmCSBV-SDLY-2016 from A. mellifera was more similar to the CSBV strains isolated than AmSBV-UK-2000. The CSBV strains isolated formed a relatively independent separation group; there was an obvious deviation compared to AmSBV-UK-2000, with maximum deviation (19.3%) occurring at 2,181 and 2,221 bp. A high degree of consistency was observed between the genomes of AmCSBV-SDLY-2016 and isolated CSBV strains, with more than 85% similarity between the complete coding regions. Among all the isolated CSBV strains, only recombinant CSBV-JLCBS-2014 was detected by RDP (Table 5) that exhibited three recombinant fragments as follows: one obtained from AmCSBV-SDLY-2016 and two derived from CSBV-FZ-2014. The minor parent was CSBV-LNQY-2008. BootScan analysis showed that the cross-recombination signals of AmCSBV-SDLY-2016 and CSBV-FZ-2014 were detected at positions 1,957–2,162 and 6,635–7,110 using AmCSBV-SDLY-2016 and CSBV-JLCBS-2014 as the query sequence, respectively (Fig. 8).