Impact of Sacbrood Virus on Larval Microbiome of Apis mellifera and Apis cerana

In this study, we examined the impact of Sacbrood virus (SBV), the cause of larval honeybee (Apis mellifera) death, producing a liquefied a larva sac, on the gut bacterial communities on two larval honeybee species, Apis mellifera and Apis cerana. SBV was added into a worker jelly food mixture and bee larvae were grafted into each of the treatment groups for 24 h before DNA/RNA extraction. Confirmation of SBV infection was achieved using quantitative reverse transcription polymerase chain reaction (RT-qPCR) and visual symptomology. The 16S rDNA was sequenced by Illumina sequencing. The results showed the larvae were infected with SBV. The gut communities of infected A. cerana larvae exhibited a dramatic change compared with A. mellifera. In A. mellifera larvae, the Illumina sequencing revealed the proportion of Gilliamella, Snodgrassella and Fructobacillus was not significantly different, whereas in A. cerana, Gilliamella was significantly decreased (from 35.54% to 2.96%), however, with significant increase in Snodgrassella and Fructobacillus. The possibility of cross-infection should be further investigated.


Introduction
The Western cavity nesting bee Apis mellifera and Asiatic cavity nesting bee Apis cerana are two domesticated bee species in worldwide apiculture [1,2]. Sacbrood virus (SBV), a common single-stranded RNA virus belonging to the genus Iflavirus in the family Iflaviridae [3], was the first virus identified in honeybees [4]. The virus infects bee larvae, causing the larval development process to fail, leading to eventual death [5]. The infected larva changes color from white to yellow and eventually dark brown, drying out after death [6]. Additionally, the virus can infect adult bees showing

Virus Inoculum
SBV was isolated from a symptomatic colony of A. mellifera larvae from the USDA honeybee research laboratory in Beltsville, Maryland. A total of 15 visually confirmed infected larvae, collected in a 50 mL tube containing 15 mL of phosphate-buffered saline (PBS), were homogenized and centrifuged at 2000× g. The supernatant was filtered using a 0.2 micron filter [27]. RNA was extracted from 200 µL of this preparation using TRIzol TM Reagent (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer's protocols. To evaluate the levels of contaminating Deformed wing virus (DWV) (the most ubiquitous virus) in our SBV inoculum, we assessed both the SBV and DWV titers. RT-qPCR confirmed high titers of SBV and relatively lower titers of DWV in the SBV-symptomatic larvae. To eliminate the contamination of DWV (and other less prevalent viruses) found in the original sample, we performed serial dilutions using PBS and retested the sample until DWV was no longer found in the sample [27]. Primers for SBV and DWV and qPCR cycling parameters can be found in reference [28].

SBV Infection in Honeybee Larvae
Fifth instar larvae were collected from three different colonies each for A. mellifera and A. cerana. Larval age was determined by size in comparison to 1-4th instars. The colonies for both species were located in Phrae, Thailand and observed to be free of SBV symptoms. Ten larvae from each of the 6 colonies were grafted using a small grafting tool and placed onto a petri dish provisioned with a worker food mixture (6% D-glucose, 6% D-fructose, 1% yeast extract and 50% royal jelly) [29] spiked with high virus titer (10 6 genome copies/larva), or no virus spike for a negative control consisting of uninoculated worker food. To assess viral infection and subsequent microbial activity, six larvae were removed from each test dish after 24 h of exposure and placed into 1.5 mL microcentrifuge with TRIzol TM Reagent for RNA and DNA extraction.

Total RNA and DNA Extraction
The larval sample stored in TRIzol TM Reagent (Thermo Fisher Scientific, Waltham, MA, USA) was extracted following the manufacturer's protocol. The two phases were extracted separately and the colorless upper aqueous phase was transferred for extraction of total RNA. The interphase of the separated phase was processed for DNA extraction using the DNeasy Blood&Tissue Kit (Qiagen, Hilden, Germany). Total RNA and DNA were quantified and qualified using a Micro-Spectrophotometer NANO-200 (Allsheng, Hangzhou Allsheng Inc., Zhejiang, China).
The phylogenetics tree of SBV was constructed by Maximum Likelihood method and present in Figures S1 and S2.

Amplicon Sequencing Analysis
The extracted DNA had at least 2 ng/µL, purity value ≥1.7 and volume ≥30 µL. The qualified DNA was sent for next-generation sequencing for 16S amplicon (V3-V4 regions) using the Illumina platform with Illumina MiSeq model by Macrogen (Macrogen Inc., Seoul, Korea). The sequences were analyzed on Miseq SOP-Mothur program [31]. The Greengenes database version 13.8 was used to compare the bacterial sequences.
Rarefaction curves were generated using OTUs and Shannon index by PAST software version 3.14 [32].

Sequences Processing
Sequences were filtered for length, 0-469 base pairs, and reduced redundancy. A redundancyreduced sequence was aligned to the Greengenes database version 13.8 as a reference and any poorly aligned sequences were removed. Chimera sequences were removed, along with any out of the ordinary sequences, before performing preliminary clustering of sequences and removing non-bacterial sequences by comparison to the Ribosomal Database Project (RDP) database.

Clustering Sequences into Operational Taxonomic Units (OTUs)
A pairwise distance between sequences was calculated, then sequences were grouped into OTUs based on the OptiClust algorithm and each OTU was assigned a name.

Alpha Diversity Analysis
Comparison of gut microbial diversities was performed using PAST version 3.14 and indexed by the Simpson index.

Beta Diversity Analysis
The non-metric multidimensional scaling (NMDS) were used to determine gut bacterial community difference between healthy and infected honeybees. To perform significance tests of beta diversity, One-way PERMANOVA was used with the Bray-Curtis similarity index.

Statistical Analysis
Differences between samples were considered significant when p < 0.05. All statistical analyses were performed using PAST version 3.14. The differences in gut communities of honeybees between the healthy and infected were carried out with an unpaired t-test with Welch's correction.
Linear discriminant analysis Effect Size (LEfSe) was used to determine the taxonomy of bacterial communities using the Galaxy application (http://huttenhower.sph.harvard.edu/galaxy) with a pairwise Wilcoxon test (<0.05) and Kruskal-Wallis test. The linear discriminant analysis (LDA) threshold score of bacteria was 2.0.
The Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) 1.0.0 (http://picrust.github.io/picrust) based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was applied to group the bacterial communities and predict the functional genes [33][34][35].

Sacbrood Virus Infection
The SBV relative expression indicating the level of infection in six bee larvae, for each treatment, was significantly higher than that in the control larvae of both A. mellifera (p = 0.001) and A. cerana (p = 0.001) ( Figure 1). In A. mellifera, the relative expression of SBV in the control and infected larvae were 3.44 ± 0.7 and 20,635 ± 3221.03, respectively. For A. cerana, the SBV expression values were 2.62 ± 0.93 in the control larvae and 7520.45 ± 1110.65 in the infected larvae. Comparison of gut microbial diversities was performed using PAST version 3.14 and indexed by the Simpson index.

Beta Diversity Analysis
The non-metric multidimensional scaling (NMDS) were used to determine gut bacterial community difference between healthy and infected honeybees. To perform significance tests of beta diversity, One-way PERMANOVA was used with the Bray-Curtis similarity index.

Statistical Analysis
Differences between samples were considered significant when p < 0.05. All statistical analyses were performed using PAST version 3.14. The differences in gut communities of honeybees between the healthy and infected were carried out with an unpaired t-test with Welch's correction.
Linear discriminant analysis Effect Size (LEfSe) was used to determine the taxonomy of bacterial communities using the Galaxy application (http://huttenhower.sph.harvard.edu/galaxy) with a pairwise Wilcoxon test (<0.05) and Kruskal-Wallis test. The linear discriminant analysis (LDA) threshold score of bacteria was 2.0.
The Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) 1.0.0 (http://picrust.github.io/picrust) based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was applied to group the bacterial communities and predict the functional genes [33][34][35].

Sacbrood Virus Infection
The SBV relative expression indicating the level of infection in six bee larvae, for each treatment, was significantly higher than that in the control larvae of both A. mellifera (p = 0.001) and A. cerana (p = 0.001) ( Figure 1). In A. mellifera, the relative expression of SBV in the control and infected larvae were 3.44 ± 0.7 and 20,635 ± 3221.03, respectively. For A. cerana, the SBV expression values were 2.62 ± 0.93 in the control larvae and 7520.45 ± 1110.65 in the infected larvae.

The Gut Bacteria of A. mellifera in the Control and Infected Larvae
The major identifiable gut bacterial component found in A. mellifera larvae was Gilliamella spp.  Tables S6 and S7).

The Gut Bacteria of A. mellifera in the Control and Infected Larvae
The major identifiable gut bacterial component found in A. mellifera larvae was Gilliamella spp.  Proteobacteria showed the highest proportion of larval microbiome in both control (84.99%) and infected larvae (91.46%) of A. mellifera followed by Firmicutes with the proportions of 11.22% in control larvae and 8.46% in infected larvae. Comparison between both honeybee species showed that there was a significant difference in proportions of Proteobacteria and Firmicutes between the control and the infected larvae ( Figure S3).
The percentage of major bacterial genera in control larvae of A. mellifera consisted of Gilliamella spp.  Tables S8 and S9).
The larval microbiome percentage in the control larvae of A. cerana consisted of Proteobacteria (80.96%) followed by Firmicutes (15.34%). For the infected, we found that the Proteobacteria (17.76%) was the highest proportion in the infected larvae followed by Firmicutes (76.15%). There was a significant difference in Proteobacteria and Firmicutes proportions between the control and the infected larvae.
The proportion of major bacterial genera in the control larvae of A. cerana was Gilliamella spp.  Proteobacteria showed the highest proportion of larval microbiome in both control (84.99%) and infected larvae (91.46%) of A. mellifera followed by Firmicutes with the proportions of 11.22% in control larvae and 8.46% in infected larvae. Comparison between both honeybee species showed that there was a significant difference in proportions of Proteobacteria and Firmicutes between the control and the infected larvae ( Figure S3).
The percentage of major bacterial genera in control larvae of A. mellifera consisted of Gilliamella spp. The larval microbiome percentage in the control larvae of A. cerana consisted of Proteobacteria (80.96%) followed by Firmicutes (15.34%). For the infected, we found that the Proteobacteria (17.76%) was the highest proportion in the infected larvae followed by Firmicutes (76.15%). There was a significant difference in Proteobacteria and Firmicutes proportions between the control and the infected larvae.
The proportion of major bacterial genera in the control larvae of A. cerana was Gilliamella spp.

Alpha Diversity Indices of Gut Bacterial Communities
Alpha diversity was calculated with Simpson parameters (Figure 3) using PAST version 3.14 (Shannon and Chao 1 parameters were showed in Figure S4). For the larval microbiome abundance between the control and the infected groups in A. mellifera, the p-values of the diversity indexed by the Simpson index was 0.192 (Figure 3a). In A. cerana, the p-values of the diversity indexed by the Simpson index was 0.028 (Figure 3b). The rarefaction plot between OTU and the Shannon index is shown in Figure S5.

Alpha Diversity Indices of Gut Bacterial Communities
Alpha diversity was calculated with Simpson parameters (Figure 3) using PAST version 3.14 (Shannon and Chao 1 parameters were showed in Figure S4). For the larval microbiome abundance between the control and the infected groups in A. mellifera, the p-values of the diversity indexed by the Simpson index was 0.192 (Figure 3a). In A. cerana, the p-values of the diversity indexed by the Simpson index was 0.028 (Figure 3b). The rarefaction plot between OTU and the Shannon index is shown in Figure S5.

Beta Diversity Analysis of Larval Bacterial Communities
The non-metric multidimensional scaling (NMDS) based on Bray-Curtis indices shows the composition of the gut bacterial communities of honeybees. The NMDS (Figure 4a,b) of A. mellifera and A. cerana were able to distinguish the gut bacteria flora between the control and SBV-infected larvae. The NMDS stress values of A. mellifera and A. cerana were 0.03 and 0.063, respectively.
The bacterial communities in infected groups in A. cerana were identified by LEfSe analysis and the differences were analyzed by one-way ANOVA. We found that unclassified Clostridiales (LDA = 2.55, p = 0.04) was affected by SBV infection, while the unclassified Orbaceae (LDA = 2.34, p = 5.53 × 10 −6 ) and Gilliamella (LDA = 2.71, p = 0.0012) were abundantly present in the control group (Figure 4c; Table S10).

Beta Diversity Analysis of Larval Bacterial Communities
The non-metric multidimensional scaling (NMDS) based on Bray-Curtis indices shows the composition of the gut bacterial communities of honeybees. The NMDS (Figure 4a,b) of A. mellifera and A. cerana were able to distinguish the gut bacteria flora between the control and SBV-infected larvae. The NMDS stress values of A. mellifera and A. cerana were 0.03 and 0.063, respectively.
The bacterial communities in infected groups in A. cerana were identified by LEfSe analysis and the differences were analyzed by one-way ANOVA. We found that unclassified Clostridiales (LDA = 2.55, p = 0.04) was affected by SBV infection, while the unclassified Orbaceae (LDA = 2.34, p = 5.53 × 10 −6 ) and Gilliamella (LDA = 2.71, p = 0.0012) were abundantly present in the control group (Figure 4c; Table S10).

Functional Gene Prediction from 16S rRNA
The groups showing the most pronounced different functional genes between the control and SBVinfected bee larvae are shown on the heat map (Figure 5a). The functional groups' biosyntheses of ansamycins (p = 6.71 × 10 −10 ) and bacterial chemotaxis (p = 6.89 × 10 −12 ) were found in A. cerana (Figure 5b). The LDA score of these functions was more than 2.0 (Figure 5c (Table  S11).

Functional Gene Prediction from 16S rRNA
The groups showing the most pronounced different functional genes between the control and SBV-infected bee larvae are shown on the heat map (Figure 5a). The functional groups' biosyntheses of ansamycins (p = 6.71 × 10 −10 ) and bacterial chemotaxis (p = 6.89 × 10 −12 ) were found in A. cerana ( Figure 5). The LDA score of these functions was more than 2.0 (Figure 5c (Table S11).

Discussion
In this study, we investigated the interaction between viral infection and the larval microbiome using an NGS technique to compare the difference in gut microbial communities between control larvae and SBV-infected larvae, in two different species of honeybee, A. mellifera and A. cerana.
There are five species of core gut bacteria that dominated in this study in honeybees (Apis sp.) including Gilliamella apicola, Snodgrassella alvi [36], Lactobacillus Firm-4, Lactobacillus Firm-5 clade [37,38] and Bifidobacterium asteroids [39,40]. These bacteria are often found in adult worker bees [41]. Lactobacillus Firm-4 and Lactobacillus Firm-5 are the two most abundant species, followed by G. apicola, S. alvi and B. asteroids. Even though B. asteroids is less common in larvae, it is commonly found in adult worker bees [41].
In this study, Gilliamella spp. was the most abundant larval microbiome, and Snodgrassella spp. was found in smaller proportions, in both A. mellifera and A. cerana. It is possible that the use of bee larvae in this study was not equivocal to what occurs in adult bees. However, for the study of this viral brood disease, it is necessary. Snodgrassella spp., which is considered to be one of two major species in the honeybee microbiome, had a lower proportion of Snodgrassella in this study than previous studies. Lactobacillus detected in this study could only be grouped to the Lactobacillus spp. and unclassified Lactobacillus, as we were unable to identify the species. Bifidobacterium is found in small amounts in adult bees, therefore, the use of bee larvae in this study may account for the absence of Bifidobacterium in our results.
Cox-Foster et al. (2007) [42] report on a metagenomic survey indicating a high relative abundance of γ-proteobacteria in hives exhibiting symptoms of colony collapse disorder (CCD) and the disappearance phenomenon of worker bees from a colony, more so than in apparently healthy hives. This research reported the presence of Firmicutes and α-Proteobacteria (related to the genus Lactobacillus and acetic acid bacteria, respectively) was dramatically reduced in infected bees. In our study, the comparison between control and infected larvae of Apis spp., agrees with Cox-Foster et al. (2007) [42] in that A. mellifera larvae shows increased detection of γ-Proteobacteria taxa (Gilliamella and Frischella) and decreasing detection of Firmicutes (Lactobacillus). On the other hand, A. cerana larvae demonstrated an increase in Firmicutes (Lactobacillus and Fructobacillus), and a decrease in α-Proteobacteria (Bombella) and γ-Proteobacteria (Gilliamella), while Frischella, which belong to γ-Proteobacteria, still showed more abundance in the SBV-infected larvae.
We demonstrated dramatic changes in the larval microbiome of A. cerana, more so than in A. mellifera when infected with SBV. It is important to note that A. mellifera is the source of the SBV used in the inoculum for study. These results show that SBV isolated from A. mellifera is infective to A. cerana as well. Reddy et al. 2016 [43] described the genome of different SBV strains compared to Korean-SBV, which shares more than 90 percent similarity. A. mellifera may have experienced host adaptation to this virus, which could explain the less severe changes found in the larval microbiome compared to those seen in A. cerana.
In addition, research by Engel et al. (2012) [44] on genomic analysis of gut microbiota revealed some members of the gut microbiome in adult bees promoted increased quality of life of the host. This research showed high abundance of pectin-degrading enzyme genes in the gut metagenome of the adult honeybee. The extreme decrease in Gilliamella in this study in A. cerana larvae may reveal that the activity in pectin degradation comes from this bacterium, supporting the theory that Gilliamella helps decompose pectin for honeybees [44]. This could negatively impact A. cerana larval longevity during infection with this virus.
The functional gene analysis base on the KEGG database allowed for an analysis of the biological activities found in the bacterial communities of each bee larval sample group. The most abundant functional gene found in A. cerana larvae was ansamycins, the group of antibiotics which have a naphthalene ring. The antibiotics belonging to the ansamycins group are rifamycins, halomycins, tolypomycin Y, streptovaricins, naphthomycin [45], and geldanamycin [46]. Rifamycins, halomycins, tolypomycin Y and streptovaricins have biological activity against Gram-positive bacteria and Mycobacterium tuberculosis by inhibiting RNA synthesis [44]. Naphthomycin has reported cytotoxicity against A-549 and P388 cell lines, but inhibitory activity against Gram-positive bacteria and M. tuberculosis was not observed [47]. The ansamycins group mostly produced by actinomycetes, Streptomyces sp., has not yet been reported in the family of Clostridiales in the phylum Firmicutes. In this analysis, biosynthesis of ansamycins function was the most abundant in infected A. cerana larvae, which relates to increasing unclassified Clostridiales. This relation may explain that Clostridiales could be able to produce ansamycin antibiotics, but under which conditions has yet to be studied. However, actinomycetes was found in a small amount in infected A. cerana larvae, which did not affect the gut communities as determined by LEfSe analysis. Therefore, we suggest that actinomycetes are not related to the functional genes. Another of the most abundant functional genes was bacterial chemotaxis, the movement toward or away from chemicals [48], that could be found in many bacteria including Clostridiales. This could explain the increased movement of bacteria after viral infection in A. cerana larvae. When bee larva become infected, this could affect the activities of organs and change the chemical composition in larvae. Subsequently, gut bacteria may be attracted to the new chemical composition or alternatively, move away from toxic chemicals produced in infected larvae, leading to increase the bacterial chemotaxis activities in infected A. cerana. The advantage and/or disadvantage of this chemotaxis against Sacbrood virus needs to be further analyzed. Nevertheless, the 16s rRNA gene prediction in this study was analyzed from adult honeybees, of which Gilliamella (and Orbaceae) had an overwhelming versatility of gene function. We suggest that information on the larval microbiome, gained in future studies, might be helpful to explain their exact function.

Conclusions
The A. mellifera-SBV used in this study led to a cross infection in A. cerana and caused significant changes in the larval microbiome, while in A. mellifera, less changes in the bacterial community were observed. It is possible that A. mellifera may have a host adaptation to this virus, or A. mellifera-SBV does not affect the bacterial community. On the contrary, in A. cerana, A. mellifera-SBV infection may or may not cause symptoms of infection, but the A. mellifera-SBV dramatically changed the larval microbiome in A. cerana. Therefore, it is imperative to separate A. cerana, in apiculture, from A. mellifera or other bees that pose the risk of being a Sacbrood virus carrier, passing a virus to an un-adapted host.
Supplementary Materials: The following are available online at http://www.mdpi.com/2075-4450/11/7/439/s1, Figure S1: Phylogenetic tree of Sacbrood virus in infected honey bees. Figure S2: Bootstrap consensus tree of Sacbrood virus in infected honey bees. Figure S3: Proportion of gut bacterial community (phylum) found in healthy and infected honeybees. Figure S4: Alpha diversity boxplots compare the gut microbiota between the control and infected worker larvae of A. mellifera (a) and A. cerana (b) designated by Shannon and Chao 1 index. Figure S5: Rarefaction curves of bacterial OTUs and Shannon index constructed by PAST software version 3.14. Table S1. Sacbrood virus sequence alignment by Basic Local Alignment Search Tool (BLAST-NCBI), to confirm SBV, extracted from 3 sampled infected larvae of each A. mellifera and A. cerana. Table S2. Beta actin sequence alignment by Basic Local Alignment Search Tool (BLAST-NCBI), to confirm the beta actin gene of A. cerana. Table S3. Beta actin sequence alignment by Basic Local Alignment Search Tool (BLAST-NCBI), to confirm the beta actin gene of A. mellifera. Table S4. Ribosomal Protein S5 sequence alignment by Basic Local Alignment Search Tool (BLAST-NCBI), to confirm the ribosomal protein S5 gene of A. cerana. Table S5. Ribosomal protein S5 sequence alignment by Basic Local Alignment Search Tool (BLAST-NCBI), to confirm the ribosomal protein S5 gene of A. mellifera. Table S6: Proportion of bacterial genera found in control A. mellifera. Table S7: Proportion of bacterial genera found in SBV infected A. mellifera. Table S8: Proportion of bacterial genera found in control A. cerana. Table S9: Proportion of bacterial genera found in SBV infected A. cerana. Table S10. LDA score of Functional gene prediction based on KEGG database. Table S11. LDA score of bacteria affected in A. cerana.