No effect of Bt Cry1Ie toxin on bacterial diversity in the midgut of the Chinese honey bees, Apis cerana cerana (Hymenoptera, Apidae)

Cry1Ie protein derived from Bacillus thuringiensis (Bt) has been proposed as a promising candidate for the development of a new Bt-maize variety to control maize pests in China. We studied the response of the midgut bacterial community of Apis cerana cerana to Cry1Ie toxin under laboratory conditions. Newly emerged bees were fed one of the following treatments for 15 and 30 days: three concentrations of Cry1Ie toxin (20 ng/mL, 200 ng/mL, and 20 μg/mL) in sugar syrup, pure sugar syrup as a negative control and 48 ng/mL imidacloprid as a positive control. The relative abundance of 16S rRNA genes was measured by Quantitative Polymerase Chain Reaction and no apparent differences were found among treatments for any of these counts at any time point. Furthermore, the midgut bacterial structure and compositions were determined using high-throughput sequencing targeting the V3-V4 regions of the 16S rDNA. All core honey bee intestinal bacterial genera such as Lactobacillus, Bifidobacterium, Snodgrassella, and Gilliamella were detected, and no significant changes were found in the species diversity and richness for any bacterial taxa among treatments at different time points. These results suggest that Cry1Ie toxin may not affect gut bacterial communities of Chinese honey bees.

Scientific RepoRts | 7:41688 | DOI: 10.1038/srep41688 Cry1Ie is a novel toxin with a different insecticidal mechanism from that of other Bt toxins currently applied for pest control purposes. This may help to delay the evolution of insect resistance to Bt toxins 29,30 . Some investigators have confirmed that Cry1Ie may be a good candidate to be produced in new commercial Bt variety for pest control [31][32][33] . However, except for two recent studies that have assessed the effect of Cry1Ie toxin on Western honey bee (Apis mellifera ligustica) under laboratory conditions 18,28 , the impacts of Cry1Ie toxin on honey bees are poorly understood.
The Asian honey bee, Apis cerana, is a honey bee species indigenous to Asia and is the second most populous honey bee species in China 34 . Compared to western honey bees (A. mellifera) which are the most abundant bee species managed in China, A. cerana has many unique characteristics that make them useful pollinators; including better adaptability to local climate, disease resistance and the ability of utilizing scattered nectar 35,36 . Thus, proper risk assessment of any new Bt toxin on this species is indispensable, especially for China. However, studies evaluating the effects of Cry1Ie toxin on Chinese honey bees, Apis cerana cerana, are scarce. Here, we used midgut bacterial community as a parameter to identify the impacts of Cry1Ie toxin on Chinese honey bees, Apis cerana cerana, through next-generation sequencing technology on the MiSeq platform.

Results
Changes in bacterial 16S rRNA gene abundance. The bacterial abundance was characterized by calculating 16S rRNA gene copy numbers using an absolute quantification qPCR method. For quantifying the copy numbers, a standard curve of y = − 0.3628x + 10.3199 (y = the logarithm of plasmid copy number to base 2, x = Ct value, R 2 = 0.995) was established as shown in Fig. S2. The copy number of 16S rRNA gene in each sample was calculated based on the standard curve (Table S1), then log-transformed prior to statistical analysis to satisfy the normality assumptions. The log10-transformed the copy numbers across different treatments are shown in Fig. 1, and these values were not significantly different among different treatments at any sample time (P > 0.05, Table S2). These results suggest that dietary Cry1Ie toxin has no significant impact on honey bee midgut bacterial abundance.
Intestinal bacterial communities of adult Chinese honey bees. In the current study, we characterized the Chinese honey bee midgut bacterial community via 16S rRNA amplicons sequencing on Illumina MiSeq platform. Paired-end sequencing of 16S rRNA V3-V4 gene produced a total of 2,071,950 raw sequences from 30 samples. After trimming the barcodes, primers and filtering chimeras, short and low-quality reads, 53,372 effective sequences were obtained, ranging from 29,264 to 64,739 per sample, and the average length of effective sequences reads was 419 bp (Table S3). All valid sequences were normalized for further analysis. However, rarefaction analyses showed that the number of observed species did not approach saturation (Fig. S3), even after 30,000 sequences, which means that more sequencing efforts were needed to capture more species.
Based on the average relative abundance, Proteobacteria (70%) was the most abundant phyla among the classified bacterial phyla. The classes with the highest abundance of bacteria were α -Proteobacteria (36.47%), Bacilli (25.41%), β -Proteobacteria (20.48%), and γ -Proteobacteria (14.92%). Commensalibacter of the family of Acetobacteraceae (28.65%), followed by Lactobacillus (25.31%), Snodgrassella (20.44%), were the 3 most abundant genera found.  Effects of the Bt Cry1Ie toxin on the midgut bacterial composition of Chinese honey bees. To determine the changes in the midgut bacterial communities of Chinese honey bees fed Cry1Ie toxin, a statistical analysis for the relative abundances of the 10 most abundant genera were conducted using one-way ANOVA (SPSS. 16.0). Among the 10 most abundant bacteria genera were Commensalibacter, Lactobacillus, Snodgrassella, Gilliamella, Frischella, Saccharibacter, Bartonella, Bifidobacterium, Citrobacter, and Bacteria. All were present in all treatment bee groups as major genera (>1%) (Table S4), and the histograms of these bacterial taxa are shown in Fig. 3. An analysis on the 10 most abundant genera revealed no significant differences in the relative abundances of any bacterial taxon with respect to the different treatments (Table S5, P > 0.05).
Furthermore, the difference of the midgut microbial community composition in different treatments was evaluated further using Hierarchical cluster analysis. Hierarchically clustered heat maps of all the abundant genera (relative abundance > 1%) are shown in Fig. 4. This figure shows that mostly all of the presented bacterial taxa here were clustered together corresponding to different treatments. These results demonstrate that the midgut bacterial composition across the different treatment groups was not significantly different.

Effects of the Bt Cry1Ie toxin on the midgut bacterial structure of Chinese honey bees.
To determine if Chinese honey bee midgut microbiota community structures were altered by the Cry1Ie toxin, the bacterial diversity was analyzed across different treatments. The bacterial diversity was characterized by calculating the alpha diversity parameters (Table S6). Two representatively alpha diversity parameters including observed species (sequencing depth) and the Shannon index (diversity indices), were selected for community richness comparison. Box plots of these richness estimators for individual sample groups are shown in Fig. 5. Comparison diversity indexes among groups were compared using a one-way ANOVA (SPSS. 16.0) with no significant differences in richness noticed across the five treatment groups (P > 0.05, Table S7).
For a better analysis of the relationships between gut microbiota community structures of the honey bees across the five treatments, the Principal Coordinate Analysis (PCoA) was performed based on the Unifrac metric. A 3-D plot of first principal component (PC1), second principal component (PC2), and third principal component (PC3) obtained by PCoA is shown in Fig. 6. The plot elucidated the characteristics of four botanical origins of the bacterial communities. The percentages are the percentage of total community variations explained by the components. On the PCoA plot, each symbol represents the gut microbiota of a sample. Sample points that are close together are more similar in community composition than those that are far apart. No significant differences in community structure were observed. The pattern revealed by PCoA were tested further using ANOSIM. Consistent with the PCoA plot, no significant impacts were noticed (Table S8). These data suggest that the midgut bacterial community structures in honey bees are not be influenced by the Bt Cry1Ie toxin.

Discussion
The potential risk of Bt toxin on the intestinal microbiota on non-target organisms is a major safety concern and requires careful consideration. Several studies have addressed the impact of Bt Cry toxin on honey bee intestinal microbiota; however, investigators in these studies mainly adopted conventional molecular methods, such as terminal restriction fragment length polymorphism (T-RFLP) 23,27 and polymerase chain reaction-denaturing gradient gel electrophoresis (PCR-DGGE) 25,26 to investigate the impacts. Many studies claimed that these conventional molecular methods were considered specific, but low-throughput and low sensitivity in characterizing microbial ecology merely provides preliminary information. Thus, the use of these methods to investigate the changes in microbial community structure and diversity cannot provide comprehensive insights to their differences and may result in some bacteria remaining undetected. For instance, previous studies have reported that Lactobacillus plays a significant role in honey bee health and that this genus is one of the core microflora within honey bees 37,38 . However, when PCR-DGGE was used, Lactobacillus was not detected during a study on the influence of transgenic cry1Ah maize pollen on the midgut bacteria community of the Chinese honey bee 25 .
As molecular biological technology continues to develop, some burgeoning technologies such as high-throughput sequencing is being applied frequently to characterize microbial ecology, including in the honey bee intestinal microbiota 20,28,[39][40][41][42] . These technologies have many benefits including high sensitivity, high accuracy and short processing times 43,44 . Moreover, these approaches have enough sequencing depth to cover complex bacterial communities, and the diversity in microbial populations is significantly higher than previously estimated by earlier conventional molecular methods 45,46 . In previous work, this method has been first applied into the Cry1Ie toxin risk assessment of another bee species A.mellifera 28 . In the current study, a 16S rRNA amplicons sequencing method was also used to examine the potential effects of Cry1Ie toxin on the midgut bacterial communities of Chinese honey bees. Our findings correspond with our previous reports that neither the midgut bacterial diversity nor their compositions were affected when the bees were exposed to Cry1Ie toxin 28 . These results suggest that Cry1Ie does not impact the gut bacterial communities of honey bees.
Honey bee intestinal microbiota has been extensively investigated using both culture-dependent and culture-independent methods. In this study, we also characterized the dominant midgut bacteria of Chinese honey bees through Illumina MiSeq paired-end sequencing of 16S rRNA gene, which were Commensalibacter, Lactobacillus, Snodgrassella, Gilliamella, Frischella, Saccharibacter, Bartonella, Bifidobacterium and Citrobacter from α -, β -, γ -Proteobacteria, Firmicutes and Actinobacteria. In agreement with previously published results 20,42 , all these conserved bacterial groups have been found in Chinese honey bees in our study. Reports indicate that adult honey bee worker guts are dominated by a few specific gut bacterial groups; merely the proportions of these phylotypes differ depending on diet and species 47,48 . For instance, a previous study investigated the intestinal microbiota of another bee species (A.mellifera), γ -Proteobacteria was the most common group of bacteria 28 ; while, in this study, the most abundant group of bacteria was α -Proteobacteria.
so the high levels of Commensalibacter in our study is probably related to the sugar diets fed to honey bees in the laboratory and may not accurately reflect the Commensalibacter levels found naturally in A. cerana.
In the current study, we used qPCR and high-throughput sequencing to assess the effects of Cry1Ie toxin on the midgut bacterial community of Chinese honey bees. No significant differences in the midgut bacterial abundance and composition were observed among the treatments and thus, we conclude that Cry1Ie does not affect bacterial diversity in the midguts of Chinese honey bees. Our work establishes a foundation leading to increased understanding of the effects of this toxin on honey bees and provides laboratory-based evidence for its future use  in transgenic plants to control insect pests. However, laboratory research is merely the first step of the biosafety assessment of a transgenic crop 53,54 . In order to provide a comprehensive risk assessment, future work under more realistic conditions (i.e., bees fed on GM pollen or field trials) is needed.

Methods
Honey bees and Cry1Ie toxin. Honey bee queens were caged on a brood comb to achieve a uniform age.
The combs containing nine-day old capped Chinese honey bee pupae, were brought from apparently healthy colonies located at the Department of Bee Protection and Biological Safety, Institute of Apicultural Research, Chinese Academy of Agricultural Sciences (CAAS), Beijing, China (40°00′ 28″ N, 116°12′ 18″ E). The brood combs were placed in the incubator at 34 ± 1 °C. All experimental work was conducted in June of 2015.
The Cry1Ie toxin used in the experiments was kindly provided by Prof. Jie Zhang from the State Key Laboratory for Biology of Plant Diseases and Insect Pests, the Institute of Plant Protection, CAAS, Beijing, China. The toxin was stored at − 20 °C and mixed thoroughly with sugar syrup (60% w/v sucrose solution) to obtain the desired concentrations in our study (20 ng/mL, 200 ng/mL, and 20 μ g/mL).
Experimental design and laboratory feeding of honey bees. Newly emerged adult bees (<6 h old) were collected from these frames and randomly placed into wooden, mesh-sided cages (dimensions of 10 cm × 7 cm × 8 cm) (Fig. S1). Each cage contained 30 bees and all cages were randomly assigned to the different treatment groups and maintained in a dark incubator (30 ± 1 °C, 60 ± 10% relative humidity) until the end of the experiment. In the study, honey bees were exposed to Cry1Ie toxin via an artificial sucrose diet. Altogether, there were five defined diets which included pure sugar syrup (60% w/v sucrose solution, the negative control), three concentrations of Cry1Ie toxin (20 ng/mL, 200 ng/mL, and 20 μ g/mL) 18 and 48 ng/mL imidacloprid syrup, the positive control 55 . There were three replicates for each treatment. A surplus of pollen and treatment syrup were fed ad libitum and were replaced daily until further honey bee dissection and DNA extraction.
Honey bee dissection and DNA extraction. After 15 and 30 d of exposure to treatments, 5 bees were randomly removed from each cage. These bees were placed in the freezer at − 20 °C for 10 s to render them inactive, and their midguts were isolated on ice using sterile forceps. Dissected midguts from the same treatment group were pooled in a 2 mL eppendorf tube and immediately frozen in liquid nitrogen for subsequent DNA extraction.
The MiSeq sequencing of 16S rRNA gene amplicons. The V3-V4 region of the bacterial 16S rRNA gene was targeted with the barcoded primer pair 341 f/806r (341 F: CCTAYGGGRBGCASCAG, 806 R: GGACTACNNGGGTATCTAAT) for the microbial community diversity analysis 57,58 . All PCR reactions in a 30 μ L mixture contained 15 μ L 2× Phusion Master Mix (New England Biolabs, USA), 1 μ L of forward and reverse primers, and about 10 ng DNA template. Amplification conditions were as follows: initial denaturation at 98 °C for 1 min, followed by 35 cycles of denaturation at 98 °C for 10 s, annealing at 50 °C for 30 s, extension at 72 °C for 30 s, and a final extension step at 72 °C for 5 min. The PCR products were checked by electrophoresis on 2% agarose gels and purified with GeneJET Gel Extraction Kit (Thermo Scientific). Finally, according to protocols described by Caporaso 58 , the purified products were sequenced on the Illumina MiSeq 250 platform (Illumina, San Diego, CA, USA) at Novogene Bioinformatics Technology Co., Ltd, Beijing, China. Data analysis. Paired-end reads were assigned to samples according to the unique barcode of each sample.
After cutting off the barcode and primer sequence, the reads that were shorter than twice the length of reads were merged into single, longer sequences using FLASH v1.2.7 (http://ccb.jhu.edu/software/FLASH/) 59 . The splicing sequences were called "raw reads". The low-quality sequences (shorter than 200 bp, average quality value of < 25) and the chimera sequences then were filtered out from downstream analysis using QIIME V1.7.0 software package (http://qiime.org/index.html) 60 and UCHIME algorithm (http://www.drive5.com/usearch/manual/uchime_ algo.html) 61 with the default parameters. The obtained sequences were called "clean reads" and normalized to make the samples compared at the same sequencing depth for the following analysis.
These normalized sequences were classified into operational taxonomic units (OTUs) at 97% similarity using UPARSE pipeline v7.0.1001 (http://drive5.com/uparse/) 62 . The taxonomy of the OTUs was assigned by blasting against SILVA SSU database 119 63,64 with default parameters. Alpha diversity included Shannon index and observed species of each sample and beta diversity included both unweighted and weighted Unifrac distances between samples, were performed with QIIME (Version 1.7.0) and displayed with R software (Version 2.15.3).
Scientific RepoRts | 7:41688 | DOI: 10.1038/srep41688 Statistical analysis. To determine whether the composition and structure of midgut bacterial communities differed significantly among treatments, statistical comparisons were made across different treatments on the 16S rRNA gene copy numbers, the dominant midgut bacterial genera composition and the richness estimators. The bacterial counts and relative abundance values were normalized with log-transformations to the base 10 prior to statistical analysis. All statistical analyses were conducted between each treatment with one-way analysis of variance (ANOVA) followed by Tukey's HSD using the software package SPSS 16.0 (IBM Co., Armonk, NY, USA).