Original Research ARTICLE
Transcriptome Analysis of Newly Emerged Honeybees Exposure to Sublethal Carbendazim During Larval Stage
- 1College of Animal Science and Technology, Yangzhou University, Yangzhou, China
- 2College of Horticulture and Plant Protection, Yangzhou University, Yangzhou, China
- 3Jiangsu Agri-animal Husbandry Vocational College, Taizhou, China
- 4Mudanjiang Branch of Heilongjiang Academy of Agricultural Sciences, Harbin, China
There are increasing concerns regarding the impact of agrochemical pesticides on non-target organisms. Pesticides could cause honeybee abnormal development in response to neurotoxins such as neonicotinoid. However, knowledge of carbendazim, a widespread fungicide in beekeeping practice, influencing on honeybee (Apis mellifera L.) brain development is lacking. Large-scale transcriptome approaches were applied to determine the changes in global gene expression in the brains of newly emerged honeybees after carbendazim exposure during the larval stage. To further understand the effects of carbendazim on the brain development of honeybees, the functions of differentially expressed genes were compared between the treatment and control groups. We found that neuroregulatory genes were down-regulated after carbendazim exposure, which suggest the neurotoxic effects of this fungicide on honeybee nervous system. Carbendazim exposure also altered the expression of genes implicated in metabolism, transport, sensor, and hormone. Notably, larvae in the carbendazim-treated group observed longer time to shift into the dormant pupal state than the control group. Moreover, a low juvenile hormone and high ecdysone titers were found in the treatment group compared to control group. The data is the first report of neurotoxic effects on honeybee caused by carbendazim, and the sublethal carbendazim may disturb honeybee development and is a potential chemical threating the honeybee colonies.
Bees are pollinators of great social, economic, and ecological importance in agricultural production and biodiversity maintenance (Potts et al., 2010). Pesticide application has considerably reduced bee colonies in many countries, and this issue is perceived as a threat to biosecurity and environmental health (Goulson et al., 2015). The use of low-toxicity pesticides has alleviated the acute effects of highly toxic pesticides to some extent, but there are increasing concerns regarding the sublethal effects of low-toxicity pesticides on non-target organisms (Desneux et al., 2007; Doublet et al., 2015; Fauser et al., 2017).
Traditionally, the classical laboratory method for estimating the side effects of chemicals on beneficial arthropods was to determine a median lethal dose (LD50) or lethal concentration. However, the estimated lethal dose during acute toxicity tests may only be a partial measure of the deleterious effects. In addition to direct mortality induced by pesticides, their sublethal effects on arthropod physiology and behavior must be considered for a complete analysis of their impact (Desneux et al., 2007). Previous studies have shown that pesticides can affect honeybee (Apis mellifera L.) behavior, neurophysiology, biochemistry, longevity, development, fecundity, immunity, orientation, mobility, learning, oviposition, and other characteristics (Papaefthimiou and Theophilidis, 2001; Desneux et al., 2007; Fischer et al., 2014; Sandrock et al., 2014; Kessler et al., 2015; Sánchez-Bayo et al., 2016; Wu et al., 2017). Different types of pesticides have distinct mechanisms and target organisms. Neurotoxicant and neonicotinoids are known to severely disturb the honeybee nervous system, including synaptic transmission, navigation, and learning ability (Jin et al., 2015; Stanley et al., 2016; Wu et al., 2017). These behavioral changes usually result in colony failure and they have raised public awareness regarding honeybees’ health and welfare. A 2-year ban on the usage of the three most common types of neonicotinoid was launched in 2013. Recently, Europe strengthened the ban on bee-harming pesticides in response to the concern that the toxicity of these chemicals actually increases over time.
The present study investigated the toxic effects of the broad-spectrum fungicide carbendazim on honeybees. In China, carbendazim has been used to control fungal pathogens on camellia, rape nectar plants, and other crops. In mainland China, the detection rate of carbendazim in 282 apple samples was 81.9%, which is one of the highest (Ye et al., 2016). Moreover, Zhou measured pesticide residues of 48 pollen samples collected from eight provinces in mainland China. As a fungicide with the highest detection rate (77.1%), the maximum concentration of carbendazim reached 4,516 ng/g.(Zhou et al., 2016). A survey conducted by the China Institute of Apicultural Research also showed that the detection rate of carbendazim was the highest among all pesticides identified in bee products. In fact, carbendazim residues have been found in pollen from many countries other than China (Mullin et al., 2010; Kasiotis et al., 2014; Sanchez-Bayo and Goka, 2014; David et al., 2015; Calatayud-Vernich et al., 2018). At present, migratory beekeeping has become a major way of commercial apiculture in China. Therefore, carbendazim exposure might be inevitable, and it might be transported from one polluted farm to others as bees forage and pollinate. However, there are few studies on the effects of carbendazim on honeybee. Hence, it is necessary to investigate the sublethal effects of carbendazim on the growth and development of honeybees.
In recent years, constant news on the collapse of bee colonies stimulated our attention (Cox-Foster et al., 2007; Luo et al., 2008; Van der Zee et al., 2012). Completion of honeybee genome sequencing in 2006 revealed that honeybees have relatively few genes encoding detoxifying enzymes (Honeybee Genome et al., 2006). Because carbendazim is frequently applied as a crop fungicide, it might be considered as a possible factor contributing to the mass death of bee colonies in China. Here, high-throughput RNA-sequencing (RNA-Seq) was utilized to identify and quantify the expression levels of the genes transcribed in the heads of adult honeybees whose larvae had been exposed to carbendazim. The resulting transcriptome library may help elucidate the sublethal effects of carbendazim on honeybees.
Materials and Methods
Honeybee and Chemicals
Honeybees (A. mellifera) were obtained from the Honeybee Research Institute of Yangzhou University, Yangzhou, China. Colonies were maintained in accordance with standard practices. Each colony had one young egg-laying queen and a population of nine comb frames with pollen, pupae, larvae, and honey. Carbendazim is distributed as a powder and it is soluble in organic solvents like dimethyl sulfoxide (DMSO). Here, DMSO was used because its amount was unimportant in the final test solution (0.1 %, v/v), and it is a widely accepted co-solvent in toxicology experimental research involving honeybees (Suchail et al., 2000; Faucon et al., 2005; Yang et al., 2012; Wu et al., 2017).
The comprehensive and effective stress concentrations of carbendazim applied here were based on published data from the survey of pesticide residues in 48 pollen samples collected from eight provinces in mainland China (Zhou et al., 2016). In this report, the maximum carbendazim concentration detected was 4.516 ng mg−1, and, therefore the rounding value of 5 ng mg−1 was adopted in the present study. Generally, different stress methods and times may result in discrepant conclusions, even when experimental subjects are exposed to identical concentrations in toxicology tests. Thus, a necessary pre-experiment was conducted to demonstrate that there was no apparent difference in mortality between control and 5 ng mg−1 carbendazim-treated groups (χ2 = 5.642, adjusted P = 0.078 > 0.05) (Supplementary Table S1), even at high concentrations. Hence, 5 ng mg−1 was considered as a reasonable and practical sublethal concentration of carbendazim to apply to honeybees in our experiment.
Sample Preparation and in vitro Artificial Larval-Rearing
The sublethal effects of carbendazim on honeybees were assessed using an artificial larva-rearing system modified from previous studies (Crailsheim et al., 2013; Wu et al., 2017). Three honeybee colonies were selected for this experiment. An excluder was used to restrict the queen to deposit eggs in a specified empty frame in each test hive for 10 h, and the same experiments were conducted in the other two test hives simultaneously. After 72 h, 1-day-old larvae were collected from combs and maintained on a diet consisting of 0.5% yeast extract, 8% fresh honey, 36.5% ultrapure water, and 55% fresh royal jelly for 1 day. This formulation was prepared at the apiary of Yangzhou University and provided in a sterile petri dish. On the 2nd day, larvae were transferred to a 48-well plate. One larva and 20 μL feed were placed in a new 48-well and the quantity of food administered was increased each day (20 μL, 30 μL, 45 μL, and 55 μL on days 3, 4, 5, and 6, respectively) (Wu et al., 2017). Once daily for four consecutive days, second instar larvae in the 48-well plate received the artificial food adulterated with 5 ng mg−1 carbendazim dissolved in 0.1% DMSO. Larvae in the non-treated group were only fed artificial food containing 0.1% DMSO. On day 7, uric acid crystal deposition in the 48-well plate was observed and the larvae were then transferred to 24-well plates fitted with sterile filter paper. Two days later, the larvae were transferred to new 24-well plates with fresh sterile filter paper. All pupae and larvae were maintained in an incubator at constant 33°C and 70% relative humidity. The brains of five newly emerged workers were dissected off and pooled, frozen in liquid nitrogen, and prepared for RNA extraction and RNA-Seq.
RNA Isolation, Library Preparation, and Sequencing
Frozen worker honeybee head samples were homogenized in TRIzol reagent (Invitrogen, Carlsbad, CA, United States), and RNA samples were pooled for cDNA library construction. An Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States), and a Qubit fluorometer (Invitrogen, United States), were used to determine RNA quantity and quality, respectively. First-strand cDNA was synthesized from mRNA after enrichment by fragmentation with buffer and oligo (dT) magnetic beads (Invitrogen, United States). In brief, random hexamer primers were used for library preparation and amplification. A QIAquick PCR extraction kit (Qiagen, Hilden, Germany), was used to purify the double-stranded cDNA, which was then eluted with EB buffer for end-repair and addition of poly (A). The Illumina HiSeqTM 2000 (Illumina, San Diego, CA, United States), was used to sequence the cDNA library and generate 125, and 150-bp paired-end reads.
Raw Data Acquisition and Statistical Analysis
A next-generation sequencing quality control (NGS QC) toolkit was used to process raw reads (Patel and Jain, 2012), removing low-quality reads and poly-N. Using Bowtie2 (Langmead and Salzberg, 2012), the resulting clean reads were mapped to the reference A. mellifera genome1 (Langmead and Salzberg, 2012). Despite sequencing discrepancy and gene length, gene expression data were normalized to reads per kilobase of transcriptome per million mapped reads (RPKM) in order to evaluate the gene expression levels (Mortazavi et al., 2008). An absolute value of log2Ratio ≥ 1 (twofold change) and false discovery ratio (FDR) ≤ 0.05 were set as significance thresholds and used to filter the differentially expressed genes (DEGs) (Anders and Huber, 2012).
Functional Analyses of DEGs
A hierarchical cluster analysis was conducted on the DEGs to identify transcript expression patterns. The hypergeometric distribution model was utilized to conduct Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) enrichment analyses (Anders and Huber, 2012; Roberts and Pachter, 2013).
Validation by Quantitative Real-Time PCR (qRT-PCR)
Nine candidate differentially regulated carbendazim-sensitive genes with different biological functions were analyzed by qRT-PCR with three biological replicates. The gene-specific primers (GSPs) used are listed in Supplementary Table S2. The RNA samples for qRT-PCR validation and for sequencing were randomly and simultaneously collected. Reactions were performed with an ABI 7500 system (Applied Biosystems, Foster City, CA, United States) and SYBR Green (Vazyme, Jiangsu, China). The housekeeping gene β-actin (AB023025) was used as the internal control (Wang et al., 2018). Transcript levels were determined in relation to that of β-actin using the 2−ΔΔCt method (Langmead and Salzberg, 2012), and evaluated for significant differences between treated and non-treated groups using an independent-samples t-test in SPSS v. 16.0 (IBM Corp., Armonk, NY, United States) considering P < 0.05 as the threshold.
Juvenile Hormone (JH) and Ecdysone (Ecd) Titers in 6-Day-Old (6 days) Honeybee Larvae
Using the same artificial larval-rearing method, 15 of 6-days honeybee larvae were randomly sampled from control and carbendazim-treated groups for JH and Ecd titer measurements, respectively. Three honeybee colonies were used as biological replicates. For a better understanding of the effect of carbendazim on JH and Ecd of honeybee larvae, the high sublethal concentration of 50 ng mg−1 was also adopted.
The test samples were triturated with liquid nitrogen and 5.0 mL methyl alcohol and 3 mL trimethylpentane were added to 1.0 g of each powdered sample. This mixture was vortexed for 5 min and then centrifuged at 10,000 × g for 15 min at 4°C. The supernatant was collected, mixed with another 3 mL trimethylpentane, and centrifuged as above. The supernatant was evaporated to dryness with a gentle nitrogen flow at 40°C. The residue was reconstituted with 1.0 mL of acetonitrile/water (2:8, v/v), and this redissolved solution was filtrated through 0.22 μm micropore cellulose membrane and analyzed by ultra high-performance liquid chromatography coupled to quadrupole time-of-flight (UHPLC-Q-TOF/MS) (Sun et al., 2016). The ACQUITY UPLC system (Waters Co., Milford, MA, United States), coupled with the hybrid Q-TOF-mass spectrophotometer SYNAPT HDMS (Waters, St. Michael, United Kingdom), was used, and separation was achieved through an ACQUITY BEH RP18 column (100 mm × 2.1 mm i.d., 3.5 μm particle size) (Waters Co., United States), at the operating flow rate of 0.3 mL min−1. Mobile phase A: aqueous solution containing 0.1% formic acid; Mobile phase B: acetonitrile containing 0.1% formic acid. The system was operated in positive electrospray ionization mode (ESI+) with the following typical parameters for mass spectrometry: capillary voltage, 3.0 kV; desolvation temperature 350°C; atomizing pressure, 35 psi; nitrogen flow, 9 L min−1; scanned range, m/z = 200–2000. The independent-sample t-test in SPSS v. 16.0 (IBM Corp., Armonk, NY, United States) was adopted for evaluating significant differences between treated and non-treated groups, considering a threshold of P < 0.05.
Raw RNA-Seq Data Analysis
Transcriptome sequencing statistics for all samples are summarized in Supplementary Table S3. Approximately, 48 million clean, single-end mRNA reads were generated. Of these,> 80% were successfully aligned to the honeybee reference genome. The expression levels of the genes identified and expressed in the heads of the newly emerged workers exposed to carbendazim or non-treated were compared based on RPKM (Trapnell et al., 2010). The reported sequencing data has been approved and assigned to SRA (Sequence Read Archive) database (SRA accession number: SRP158178). Correlation analyses among the three biological replicates of each sample (Supplementary Table S4) indicated high sequencing reliability.
Screening and Functional Analyses of DEGs
Based on P < 0.05 and log2Ratio ≥ 1, 247 DEGs were identified between the carbendazim-treated and non-treated (control) samples. Of these, 187 were upregulated and 60 were downregulated in the carbendazim-treated group compared to the control group (Supplementary Table S5).
The biological functions of the 247 DEGs were classified according to the GO database, and were organized into three main categories: molecular function, biological process, and cellular component. The 692 significantly enriched GO terms successfully mapped (Q < 0.05), and the top 10 enriched terms per category were selected to further elucidate the functions of the DEGs. These included the activities of lipid transporter, monooxygenase, and oxidoreductase, as well as cuticle composition and heme-binding in the molecular function category, and neuropeptide signaling pathway, response to reactive oxygen species, fatty acid metabolism, brain development, visual perception, and drug response in the biological process category (Figure 1). The biochemical pathways of the DEGs were investigated using the KEGG database. The top 20 enriched KEGG pathways included fatty acid valine, leucine, and isoleucine degradation, tryptophan, glycine, serine, and threonine metabolism, and eukaryotic ribosome biogenesis (Figure 2).
FIGURE 1. Gene ontology (GO) top 30 terms enriched by differentially expressed genes (DEGs). The results are summarized in three main categories: biological process, cellular component, and molecular function. The x-axis indicates the second term of GO and the y-axis indicates gene percentage.
FIGURE 2. Top 20 significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of DEGs. Color and gene ratio indicate the P-value and the ratio of genes within each pathway, respectively. The x-axis indicates the enrichment score of KEGG and the y-axis indicates the name of the top 20 pathways.
Validation of Target DEGs by qRT-PCR
The expression patterns of the nine DEGs selected for validation by qRT-PCR were consistent with that from RNA-Seq (Figure 3 and Supplementary Table S2). Relative expression levels of these DEGs provided a better comparison and the highest expression levels were set to 1 (Wang et al., 2018). Their agreement indicated that the abundance of the Illumina reads closely mirrored the actual expression levels of the DEGs.
FIGURE 3. Joint analysis of RNA sequencing (RNA-Seq) and quantitative real-time PCR (qRT-PCR). The x-axis represents the functional gene name and the y-axis the relative level of their expression by RNA-Seq and qRT-PCR. Green and blue represent the values obtained from RNA-seq and qRT-PCR methods, respectively. All genes with significantly different expression between the carbendazim-treated and control groups are indicated using asterisks: ∗P < 0.05, ∗∗P < 0.01 (independent-samples t-test).
Titer of JH and Ecd in 6 Days Larvae After Carbendazim Exposure
Notably, larvae in the carbendazim-treated group needed more time to shift into the dormant pupal state than larvae in the non-treated group, which was considered abnormal we observed. The JH and Ecd titers of 6 days honeybee larvae are shown in Figure 4. As depicted, the JH titers of 6 days larvae treated with carbendazim increased by 1.36 and 1.56 times in relation to that of control larvae (Figure 4A), respectively (P < 0.05). The Ecd titers showed an opposite pattern (Figure 4B), as that of the treated group were lower than that of the control group (82%), especially in honeybee larvae exposed to the highest carbendazim stress concentration (50 ng mg−1; P < 0.05). These results indicated that the tested concentrations of carbendazim can promote JH synthesis and inhibit the secretion of Ecd in 6 days larvae.
FIGURE 4. Titers of juvenile hormone (JH; A) and ecdysone (Ecd; B) in 6-day-old larvae after carbendazim exposure and in the control group. Data presented in the table are mean ± standard error. The x-axis represents the control group [larvae fed artificial food containing 0.1% dimethyl sulfoxide (DMSO)] and the treated groups (larvae fed with artificial food adulterated with 5 or 50 ng mg−1 of carbendazim dissolved in 0.1% DMSO). Different letters above bars represent significant differences (independent-samples t-test, P < 0.05).
Rearing honeybee larvae in vitro is of great importance for research on pathogens and risk assessment (Crailsheim et al., 2013). In the present study, honeybees were exposed to the fungicide carbendazim to assess and evaluate its sublethal effects. The examined honeybees were raised in vitro from the larval stage in order to eliminate the effect of the interfering factors associated with hive rearing. Suitable carbendazim post-exposure risk evaluation samples were then prepared for RNA-Seq and qRT-PCR. Our previous in vitro larval-rearing revealed a 75 emergence rate, which was acceptable (Brodschneider et al., 2009; Aupinel et al., 2010; Crailsheim et al., 2013; Wu et al., 2017), and demonstrated the reliability of the present results.
Among the 247 DEGs identified after carbendazim exposure (5 ng mg−1), using a rigorous quantitative strategy via RNA-Seq, 187 were upregulated and 60 were downregulated in the treated group compared with the control group. The DEGs found in the present study reflect the multifaceted influence of carbendazim on honeybee physiology. Interestingly, results showed that carbendazim exposure at the larval stage significantly changed the expression of genes encoding major royal jelly proteins (MRJPs) in adult honeybee brain. In fact, MRJP1, MRJP3, MRJP4, MRJP5, and MRJP6 were upregulated in the carbendazim-exposed group (Supplementary Figure S1). A previous study showing the expression of MRJP1–9 in honeybee brains suggested that these genes might have other functions in addition to regulating royal jelly biosynthesis (Buttstedt et al., 2013). For example, MRJP1 has antimicrobial properties (Fontana et al., 2004; Rosmilah et al., 2008). To the best of our knowledge, the present study is the first to show that carbendazim exposure in honeybee upregulates MRJPs expression, a finding that merits further investigation.
In 187 upregulated genes library, 59 of these genes (approximately 31.5%), were grouped into receptor and transporter categories. Among these 59 genes, 43 were upregulated and 16 were downregulated (Supplementary Figure S1). The upregulation of these genes in response to carbendazim exposure might induce toxicant metabolism (Wu et al., 2017). For example, carbendazim stress upregulated the expression of cytochrome P450s, such as P450 6AS5 (GB49890) and 4aa1 (GB49626). The members of the CYP450 family 4 are responsible for the metabolism of endogenous substrates, including pheromones and hormones, whereas members of the CYP450 family 6 are involved in xenobiotic detoxification and transport (Corona and Robinson, 2006). On the other hand, one of the 16 transporter and receptor genes downregulated by carbendazim exposure encoded apolipoprotein D (GB54020). Therefore, this fungicide may adversely affect honeybee receptor-associated lipid metabolism. Other functional genes like G-protein coupled receptor (GB48344), which have been demonstrated to enhance spatial recognition memory in rats, was also downregulated, and this suppression might cause orientation and identification disorders in honeybee after carbendazim exposure (Riedel et al., 2003; Hawley et al., 2014). Because, both receptor and transporter proteins are essential for metabolism and development the changes in the expression of the genes encoding them might cause adverse effects, ultimately destroying the physiological system.
Transcriptome analyses indicated that sublethal carbendazim exposures in honeybees either upregulated or downregulated genes involved in immunity, detoxification, sensory processing, brain development, and metabolism. The variability of gene expression in response to carbendazim exposure suggests that this fungicide might influence multiple aspects of honeybee physiology, threatening their survival. For instance, in carbendazim-exposed honeybees the number of transcripts of Cyp4aa1-like, Cyp6as5, abaecin (GB18323), and defensin 1 (GB41428), which are involved in detoxification and immunity, was 31-, 2.8-,13-, and 197.7-fold that in control honeybees, respectively, suggesting that defense was initiated upon fungicide exposure. Importantly, certain genes associated with brain development were significantly downregulated after carbendazim exposure. These included synaptotagmin 20 (Syt20, GB54319), cyclin-dependent kinase 5 activator 1 (Cdk5alpha, GB55798), dopamine transporter (DAT, GB40867), SIF amide (GB40093), and slit homolog 1 (GB43897) (Supplementary Figure S1).
Synaptotagmins are synaptic proteins, which play important roles in calcium-mediated, depolarization-induced exocytosis, neurotransmitter release, and hormone secretion. These processes are essential for the transmission of neural signals and brain development (Fernandez-Chacon et al., 2001; Chang et al., 2018). Because synaptotagmins share a similar domain structure and a high degree of homology (Huang et al., 2017; Saheki, 2017), Syt20 might harbor similar functionality although its characteristics have not been studied to date. The specific activator of cyclin-dependent kinase 5 (CDK5), Cdk5alpha, is an important molecule for nervous system development and function and for neurodegenerative disorder pathogenesis (Kobayashi et al., 2014; Papadopoulou et al., 2015). Mutations in Cdk5 have been associated with neuronal injury and neurodegeneration and lead to the failure of synaptic cognition and modeling (Lai et al., 2015; Yousuf et al., 2016). Dopamine is an important neurotransmitter in both invertebrate and vertebrate nervous systems. It is widely distributed in the honeybee brain and it is associated with behavior regulation and olfactory memory (Schürmann et al., 1989; Beggs and Mercer, 2009; Mustard et al., 2010). Dopamine transporter removes dopamine from the synaptic cleft and deposits it into the surrounding cells. Therefore, DAT terminates the neurotransmitter signal and helps elucidate the primary mechanism of dopamine. Previous studies have shown that the queen regulates worker behavior and physiology via certain molecular mechanisms, especially dopamine signal pathways (Beggs et al., 2007; Beggs and Mercer, 2009). Therefore, DAT downregulation may affect honeybee brain development and can induce rebellious behavior in workers and colony disorder. Evidence indicates that, under axonal guidance, slit homolog 1 participates in nervous system development, especially cortex development (Hao et al., 2001; Farmer et al., 2008). The highly conserved neuropeptide SIF amide plays an important role in the olfactory neural pathway and has a neuromodulatory role in combining tactile, olfactory, and visual inputs (Verleyen et al., 2004, 2009). In the present study, the expression of the genes encoding these proteins related to brain development was validated by qRT-PCR (Figure 3 and Supplementary Table S2). The significant agreement of RNA-Seq and qRT-PCR results indicated that Illumina sequencing results closely mirrored the actual expression level of these genes and emphasized the negative effect of carbendazim on honeybee brain. In the present study, GO enrichment analysis showed that the expression of the genes involved in brain development and on neuropeptide signaling pathway was significantly decreased in honeybees exposed to carbendazim (GO:0007420 and GO:0007218). These fluctuations in gene expression may impair honeybee brain function and possibly lead to the collapse of the whole colony.
It is widely accepted that exogenous adverse factors can delay insect development, as found in honeybee (Czoppelt and Rembold, 1988; Wu et al., 2014). In the present study, we found that carbendazim exposure delays insect development. Larvae in carbendazim-treated groups needed more time to shift into the dormant pupal state and this abnormality was more obvious with increasing carbendazim concentrations (50 ng mg−1). The degree of honeybee development is often regulated by hormones, mainly JH and Ecd (Corona et al., 2007; Schulz et al., 2014). The titers of JH and Ecd in 6 days larvae found here (Figure 3) were higher and lower, respectively, in the carbendazim-treated group than in the control group, thus confirming that this fungicide accounted for the delay in development. The Kruppel-like homolog 1 (Kr-h1), a hormone regulatory gene, was one of the DEGs found in the present study. It is a key regulator of insect molting and metamorphosis and a major effector in JH signaling (Pecasse et al., 2000; Minakuchi et al., 2008). During insect development, Kr-h1 functions as a transcriptional repressor on the neurogenesis of the mushroom body and photoreceptor maturation, which is considered an interaction with the JH signaling (Shi et al., 2007; Fichelson et al., 2012; Kayukawa et al., 2012). In Drosophila melanogaster, Kr-h1 delays larval development and alters lipid metabolism (Ping et al., 2017). Furthermore, previous research reported that structural cuticle genes were downregulated as the insect entered the diapause maintenance phase of diapause development (Yocum et al., 2009). In the present study, members of a cuticle protein family: apidermin 1 (GB53115), apidermin 2 (GB53119), apidermin 3 (GB53114), cuticular protein 5 (GB40299), and cuticular protein 19 (GB50453) tended to be downregulated in the brains of newly emerged adult honeybees after carbendazim exposure (Supplementary Figure S1). This downregulation of cuticular genes might be a morphological manifestation of the delay in development. Thus, carbendazim-induced Kr-h1 upregulation may cause JH disorder and delay honeybee development.
Some other functional genes were found to change their expression patterns after carbendazim exposure. Kin of IRRE-like protein 2 (Kirrel2, GB45525) encodes pancreatic islet β-cell-expressed protein. It is an integral part of the extracellular adhesion scaffold required for normal basal insulin secretion (Yesildag et al., 2015). Insect insulin-like peptides and their signal pathways play important roles in metabolism, longevity, growth, and development (Antonova et al., 2012; Fridell, 2013; Yesildag et al., 2015). Queen brain-selective protein-1 (QBP-1, GB45741) plays a role in determining honeybee longevity (Grönke and Partridge, 2010), and its isoforms function in phase transition (Vanhems et al., 1990; Claeys et al., 2005) in locusts. Therefore, downregulation of these genes might reflect the impairment of certain unmanifested functions.
Overall, the results obtained here demonstrate that sublethal doses of carbendazim hinder growth and development and may destabilize and impede the development of honeybee colonies. Thus, they should be used to inform beekeepers and farmers about the effects of carbendazim on honeybees and provide guidance in using this fungicide carbendazim in agriculture reasonably.
This study was carried out in accordance with the recommendations of ‘Guidelines for laboratory animals of Yangzhou university, Yangzhou university animal welfare committee’ with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the ‘Yangzhou university animal welfare committe.’
KW designed the study and carried out the data analysis. R-LF and W-NJ participated in drafting the manuscript. W-WZ, X-MC, and SW participated in sample collection. LY and TJ provided advice on data analyses and helped to drafting the manuscript. F-CG and G-HC were responsible for answering questions and language Edit in the revision process.
This work was supported by grants from the National Natural Science Foundation of China (No. 31372382). Natural Science Foundation of Heilongjiang Province (No. C2017062). The Earmarked Fund for Modern Agro-Industry Technology Research System from the Ministry of Agriculture of China (CARS-45-SYZ6) to TJ.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We thank Prof. En-cheng Yang in Taiwan University, and Prof. Bao-hua Xu and Dr. Wei-xin Zhang in Shandong Agricultural University, China, for the guidance on raising honeybee in vitro, and Dr. Zheng-guo Liu in Shandong Agricultural University and Ling Pan in Shanghai Oebiotech Co., Ltd., China for advice on data analyses and technical assistance.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2018.00426/full#supplementary-material
FIGURE S1 | Correlation coefficient of samples. A heat map representation of the correlation matrix of sequencing samples.
TABLE S1 | Survival Rates. Survival rates of control and 5 ng mg−1 carbendazim-treated groups (Chi-Square (χ2 Tests).
TABLE S2 | Gene-specific primers. Specific primers of functional gene used for qRT-PCR.
TABLE S3 | Summary of mapping result of RNA-seq. All genome mapping statistics.
TABLE S4 | All differentially expressed genes. An absolute value of log2Ratio ≥ 1 (two-fold change) and false discovery ratio (FDR) ≤ 0.05 were set as significance thresholds and used to filter the differentially expressed genes (DEGs).
TABLE S5 | Functional gene in DEGs. Different kinds of functional gene among DEGs, including major royal jelly proteins, transporter and receptor, cuticle structure and Brain development related genes.
- ^ ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/002/195/GCF_000002195.4_Amel_4.5/GCF_000002195.4_Amel_4.5_genomic.fna.gz
Antonova, Y., Arik, A. J., Moore, W., Riehle, M. A., and Brown, M. R. (2012). “Insulin-like peptides: structure, signaling, and function,” in Insect Endocrinology, ed. L. I. Gilbert (New York, NY: Elsevier), 63–92.
Aupinel, P., Fortini, D., Michaud, B., Medrzycki, P., Padovani, E., Przygoda, D., et al. (2010). Honey bee brood ring-test: method for testing pesticide toxicity on honeybee brood in laboratory conditions. Julius Kühn Archiv 423,96–102.
Beggs, K. T., Glendining, K. A., Marechal, N. M., Vergoz, V., Nakamura, I., Slessor, K. N., et al. (2007). Queen pheromone modulates brain dopamine function in worker honey bees. Proc. Natl. Acad. Sci. U.S.A. 104, 2460–2464. doi: 10.1073/pnas.0608224104
Buttstedt, A., Moritz, R. F., and Erler, S. (2013). More than royal food-Major royal jelly protein genes in sexuals and workers of the honeybee Apis mellifera. Front. Zool. 10, 72–82. doi: 10.1186/1742-9994-10-72
Calatayud-Vernich, P., Calatayud, F., Simó, E., and Picó, Y. (2018). Pesticide residues in honey bees, pollen and beeswax: assessing beehive exposure. Environ. Pollut. 241, 106–114. doi: 10.1016/j.envpol.2018.05.062
Chang, S., Trimbuch, T., and Rosenmund, C. (2018). Synaptotagmin-1 drives synchronous Ca 2+-triggered fusion by C 2 B-domain-mediated synaptic-vesicle-membrane attachment. Nat. Neurosci. 21, 33–40. doi: 10.1038/s41593-017-0037-5
Claeys, I., Simonet, G., Breugelmans, B., Van Soest, S., Franssens, V., Sas, F., et al. (2005). Quantitative real-time RT-PCR analysis in desert locusts reveals phase dependent differences in neuroparsin transcript levels. Insect Mol. Biol. 14, 415–422. doi: 10.1111/j.1365-2583.2005.00572.x
Corona, M., Velarde, R. A., Remolina, S., Moran-Lauter, A., Wang, Y., Hughes, K. A., et al. (2007). Vitellogenin, juvenile hormone, insulin signaling, and queen honey bee longevity. Proc. Natl. Acad. Sci. U.S.A. 104, 7128–7133. doi: 10.1073/pnas.0701909104
Cox-Foster, D. L., Conlan, S., Holmes, E. C., Palacios, G., Evans, J. D., Moran, N. A., et al. (2007). A metagenomic survey of microbes in honey bee colony collapse disorder. Science 318, 283–287. doi: 10.1126/science.1146498
Crailsheim, K., Brodschneider, R., Aupinel, P., Behrens, D., Genersch, E., Vollmann, J., et al. (2013). Standard methods for artificial rearing of Apis mellifera larvae. J. Apic. Res. 52, 1–16. doi: 10.7717/peerj.3858
David, A., Botías, C., Abdul-Sada, A., Goulson, D., and Hill, E. M. (2015). Sensitive determination of mixtures of neonicotinoid and fungicide residues in pollen and single bumblebees using a scaled down QuEChERS method for exposure assessment. Anal. Bioanal. Chem. 407, 8151–8162. doi: 10.1007/s00216-015-8986-6
Doublet, V., Labarussias, M., Miranda, J. R., Moritz, R. F., and Paxton, R. J. (2015). Bees under stress: sublethal doses of a neonicotinoid pesticide and pathogens interact to elevate honey bee mortality across the life cycle. Environ. Microbiol. 17, 969–983. doi: 10.1111/1462-2920.12426
Farmer, W. T., Altick, A. L., Nural, H. F., Dugan, J. P., Kidd, T., Charron, F., et al. (2008). Pioneer longitudinal axons navigate using floor plate and Slit/Robo signals. Development 135, 3643–3653. doi: 10.1242/dev.023325
Faucon, J. P., Aurières, C., Drajnudel, P., Mathieu, L., Ribiere, M., Martel, A. C., et al. (2005). Experimental study on the toxicity of imidacloprid given in syrup to honey bee (Apis mellifera) colonies. Pest. Manag. Sci. 61, 111–125. doi: 10.1002/ps.957
Fauser, A., Sandrock, C., Neumann, P., and Sadd, B. M. (2017). Neonicotinoids override a parasite exposure impact on hibernation success of a key bumblebee pollinator. Ecol. Entomol. 42, 306–314. doi: 10.1111/een.12385
Fernandez-Chacon, R., Königstorfer, A., Gerber, S. H., García, J., Matos, M. F., Stevens, C. F., et al. (2001). Synaptotagmin I functions as a calcium regulator of release probability. Nature 410, 41–49. doi: 10.1038/35065004
Fichelson, P., Brigui, A., and Pichaud, F. (2012). Orthodenticle and Kruppel homolog 1 regulate Drosophila photoreceptor maturation. Proc. Natl. Acad. Sci. U.S.A. 109, 7893–7898. doi: 10.1073/pnas.1120276109
Fischer, J., Mueller, T., Spatz, A.-K., Greggers, U., Gruenewald, B., and Menzel, R. (2014). Neonicotinoids interfere with specific components of navigation in honeybees. PLoS One 9:e91364. doi: 10.1371/journal.pone.0091364
Fontana, R., Mendes, M. A., de Souza, B. M., Konno, K., César, L. M., Malaspina, O., et al. (2004). Jelleines: a family of antimicrobial peptides from the Royal Jelly of honeybees (Apis mellifera). Peptides 25, 919–928. doi: 10.1016/j.peptides.2004.03.016
Goulson, D., Nicholls, E., Botías, C., and Rotheray, E. L. (2015). Bee declines driven by combined stress from parasites, pesticides, and lack of flowers. Science 347:1255957. doi: 10.1126/science.1255957
Grönke, S., and Partridge, L. (2010). “The functions of insulin-like peptides in insects,” in IGFs: Local Repair and Survival Factors Throughout Life Span, ed. C. Yves, (New York, NY: Springer), 105–124. doi: 10.1007/978-3-642-04302-4_9
Hao, J. C., Timothy, W. Y., Fujisawa, K., Culotti, J. G., Gengyo-Ando, K., Mitani, S., et al. (2001). C. elegans slit acts in midline, dorsal-ventral, and anterior-posterior guidance via the SAX-3/Robo receptor. Neuron 32, 25–38. doi: 10.1016/S0896-6273(01)00448-2
Hawley, W. R., Grissom, E. M., Moody, N. M., Dohanich, G. P., and Vasudevan, N. (2014). Activation of G-protein-coupled receptor 30 is sufficient to enhance spatial recognition memory in ovariectomized rats. Behav. Brain Res. 262, 68–73. doi: 10.1016/j.bbr.2014.01.006
Huang, R., Zhao, J., Liu, J., Wang, Y., Han, S., and Zhao, H. (2017). Genome-wide analysis and expression profiles of NTMC2 family genes in Oryza sativa. Gene 637, 130–137. doi: 10.1016/j.gene.2017.09.046
Jin, N., Klein, S., Leimig, F., Bischoff, G., and Menzel, R. (2015). The neonicotinoid clothianidin interferes with navigation of the solitary bee Osmia cornuta in a laboratory test. J. Exp. Biol. 218, 2821–2825. doi: 10.1242/jeb.123612
Kasiotis, K. M., Anagnostopoulos, C., Anastasiadou, P., and Machera, K. (2014). Pesticide residues in honeybees, honey and bee pollen by LC–MS/MS screening: reported death incidents in honeybees. Sci. Total Environ. 485, 633–642. doi: 10.1016/j.scitotenv.2014.03.042
Kayukawa, T., Minakuchi, C., Namiki, T., Togawa, T., Yoshiyama, M., Kamimura, M., et al. (2012). Transcriptional regulation of juvenile hormone-mediated induction of Krüppel homolog 1, a repressor of insect metamorphosis. Proc. Natl. Acad. Sci. U.S.A. 109, 11729–11734. doi: 10.1073/pnas.1204951109
Kessler, S. C., Tiedeken, E. J., Simcock, K. L., Derveau, S., Mitchell, J., Softley, S., et al. (2015). Bees prefer foods containing neonicotinoid pesticides. Nature 521, 74–76. doi: 10.1038/nature14414
Kobayashi, H., Saito, T., Sato, K., Furusawa, K., Hosokawa, T., Tsutsumi, K., et al. (2014). Phosphorylation of cyclin-dependent kinase 5 (Cdk5) at Tyr-15 is inhibited by Cdk5 activators and does not contribute to the activation of Cdk5. J. Biol. Chem. 289, 19627–19636. doi: 10.1074/jbc.M113.501148
Lai, K.-O., Liang, Z., Fei, E., Huang, H., and Ip, N. Y. (2015). Cyclin-dependent kinase 5 (Cdk5)-dependent phosphorylation of p70 ribosomal S6 kinase 1 (S6K) is required for dendritic spine morphogenesis. J. Biol. Chem. 290, 14637–14646. doi: 10.1074/jbc.M114.627117
Minakuchi, C., Zhou, X., and Riddiford, L. M. (2008). Krüppel homolog 1 (Kr-h1) mediates juvenile hormone action during metamorphosis of Drosophila melanogaster. Mech. Dev. 125, 91–105. doi: 10.1016/j.mod.2007.10.002
Mullin, C. A., Frazier, M., Frazier, J. L., Ashcraft, S., Simonds, R., and Pettis, J. S. (2010). High levels of miticides and agrochemicals in North American apiaries: implications for honey bee health. PLoS One 5:e9754. doi: 10.1371/journal.pone.0009754
Mustard, J. A., Pham, P. M., and Smith, B. H. (2010). Modulation of motor behavior by dopamine and the D1-like dopamine receptor AmDOP2 in the honey bee. J. Insect Physiol. 56, 422–430. doi: 10.1016/j.jinsphys.2009.11.018
Papadopoulou, A., Siamatras, T., Delgado-Morales, R., Amin, N., Shukla, V., Zheng, Y., et al. (2015). Acute and chronic stress differentially regulate cyclin-dependent kinase 5 in mouse brain: implications to glucocorticoid actions and major depression. Transl. Psychiatry 5:e578. doi: 10.1038/tp.2015.72
Papaefthimiou, C., and Theophilidis, G. (2001). The cardiotoxic action of the pyrethroid insecticide deltamethrin, the azole fungicide prochloraz, and their synergy on the semi-isolated heart of the bee Apis mellifera macedonica. Pest. Biochem. Physiol. 69, 77–91. doi: 10.1006/pest.2000.2519
Pecasse, F., Beck, Y., Ruiz, C., and Richards, G. (2000). Kruppel-homolog, a stage-specific modulator of the prepupal ecdysone response, is essential for Drosophila metamorphosis. Dev. Biol. 221, 53–67. doi: 10.1006/dbio.2000.9687
Ping, K., Chang, K., Liu, Y., Bouska, M., Birnbaum, A., Karashchuk, G., et al. (2017). Drosophila Kruppel homolog 1 represses lipolysis through interaction with dFOXO. Sci. Rep. 7:16369. doi: 10.1038/s41598-017-16638-1
Potts, S. G., Biesmeijer, J. C., Kremen, C., Neumann, P., Schweiger, O., and Kunin, W. E. (2010). Global pollinator declines: trends, impacts and drivers. Trends Ecol. Evol. 25, 345–353. doi: 10.1016/j.tree.2010.01.007
Saheki, Y. (2017). “Endoplasmic reticulum–plasma membrane crosstalk mediated by the extended synaptotagmins,” in Organelle Contact Sites, eds Tagaya, Mitsuo, Simmen, and Thomas (New York, NY: Springer), 83–93.
Sánchez-Bayo, F., Goulson, D., Pennacchio, F., Nazzi, F., Goka, K., and Desneux, N. (2016). Are bee diseases linked to pesticides?—A brief review. Environ. Int. 89, 7–11. doi: 10.1016/j.envint.2016.01.009
Sandrock, C., Tanadini, L. G., Pettis, J. S., Biesmeijer, J. C., Potts, S. G., and Neumann, P. (2014). Sublethal neonicotinoid insecticide exposure reduces solitary bee reproductive success. Agric. For. Entomol. 16, 119–128. doi: 10.1111/afe.12041
Schulz, D. J., Pankiw, T., Fondrk, M. K., Robinson, G. E., and Page, R. E. (2014). Comparisons of juvenile hormone hemolymph and octopamine brain titers in honey bees (Hymenoptera: Apidae) selected for high and low pollen hoarding. Ann. Entomol. Soc. Am. 97, 1313–1319. doi: 10.1603/0013-8746(2004)097[1313:COJHHA]2.0.CO;2
Shi, L., Lin, S., Grinberg, Y., Beck, Y., Grozinger, C. M., Robinson, G. E., et al. (2007). Roles of Drosophila Kruppel-homolog 1 in neuronal morphogenesis. Dev. Neurobiol. 67, 1614–1626. doi: 10.1002/dneu.20537
Stanley, D. A., Russell, A. L., Morrison, S. J., Rogers, C., and Raine, N. E. (2016). Investigating the impacts of field-realistic exposure to a neonicotinoid pesticide on bumblebee foraging, homing ability and colony growth. J. Appl. Ecol. 53, 1440–1449. doi: 10.1111/1365-2664.12689
Suchail, S., Guez, D., and Belzunces, L. P. (2000). Characteristics of imidacloprid toxicity in two Apis mellifera subspecies. Environ. Toxicol. Chem. 19,1901–1905. doi: 10.1897/1551-5028(2000)019<1901:COITIT>2.3.CO;2
Sun, F., Yang, S., Zhang, H., Zhou, J., Li, Y., Zhang, J., et al. (2016). Comprehensive analysis of tiamulin metabolites in various species of farm animals using ultra-high-performance liquid chromatography coupled to quadrupole/time-of-flight. J. Agric. Food Chem. 65, 199–207. doi: 10.1021/acs.jafc.6b04377
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., Van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621
Van der Zee, R., Pisa, L., Andonov, S., Brodschneider, R., Charrière, J.-D., Chlebo, R., et al. (2012). Managed honey bee colony losses in Canada, China, Europe, Israel and Turkey, for the winters of 2008–9 and 2009–10. J. Apic. Res. 51, 100–114. doi: 10.3896/IBRA.22.214.171.124
Verleyen, P., Huybrechts, J., Baggerman, G., Van Lommel, A., De Loof, A., and Schoofs, L. (2004). SIFamide is a highly conserved neuropeptide: a comparative study in different insect species. Biochem. Biophys. Res. Commun. 320, 334–341. doi: 10.1016/j.bbrc.2004.05.173
Wang, K., Liu, Z.-G., Pang, Q., Zhang, W.-W., Chen, X.-M., Fan, R.-L., et al. (2018). Investigating the regulation of hypopharyngeal gland activity in honeybees (Apis mellifera carnica) under overwintering conditions via morphologic analysis combined with iTRAQ-Based comparative proteomics. Ann. Entomol. Soc. Am. 111, 127–135. doi: 10.1093/aesa/say012
Wu, M. C., Chang, Y. W., Lu, K. H., and Yang, E. C. (2017). Gene expression changes in honey bees induced by sublethal imidacloprid exposure during the larval stage. Insect Biochem. Mol. Biol. 88, 12–20. doi: 10.1016/j.ibmb.2017.06.016
Wu, Y. Y., Zhou, T., Wubie, A. J., Wang, Q., Dai, P. P., and Jia, H. R. (2014). Apoptosis in the nerve cells of adult honeybee Apis mellifera ligustica brain induced by imidacloprid. Acta Entomol. Sin. 57, 194–203.
Yang, E. C., Chang, H. C., Wu, W.-Y., and Chen, Y. W. (2012). Impaired olfactory associative behavior of honeybee workers due to contamination of imidacloprid in the larval stage. PLoS One 7:e49472. doi: 10.1371/journal.pone.0049472
Yesildag, B., Bock, T., Herrmanns, K., Wollscheid, B., and Stoffel, M. (2015). Kin of IRRE-like protein 2 is a phosphorylated glycoprotein that regulates basal insulin secretion. J. Biol. Chem. 290, 25891–25906. doi: 10.1074/jbc.M115.684704
Yocum, G. D., Rinehart, J. P., and Larson, M. L. (2009). Down-regulation of gene expression between the diapause initiation and maintenance phases of the Colorado potato beetle. Leptinotarsa decemlineata (Coleoptera: Chrysomelidae). Eur. J. Entomol. 106, 471–476. doi: 10.14411/eje.2009.059
Yousuf, M. A., Tan, C., Torres-Altoro, M. I., Lu, F. M., Plautz, E., and Zhang, S. (2016). Involvement of aberrant cyclin-dependent kinase 5/p25 activity in experimental traumatic brain injury. J. Neurochem. 138, 317–327. doi: 10.1111/jnc.13620
Keywords: honeybee, sublethal dose effects, carbendazim, brain development, transcriptome
Citation: Wang K, Fan R-L, Ji W-N, Zhang W-W, Chen X-M, Wang S, Yin L, Gao F-C, Chen G-H and Ji T (2018) Transcriptome Analysis of Newly Emerged Honeybees Exposure to Sublethal Carbendazim During Larval Stage. Front. Genet. 9:426. doi: 10.3389/fgene.2018.00426
Received: 04 July 2018; Accepted: 10 September 2018;
Published: 08 October 2018.
Edited by:Jianke Li, Institute of Apiculture Research (CAAS), China
Reviewed by:Chunsheng Hou, Institute of Apiculture Research (CAAS), China
Linsheng Yu, Anhui Agricultural University, China
Copyright © 2018 Wang, Fan, Ji, Zhang, Chen, Wang, Yin, Gao, Chen and Ji. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ting Ji, firstname.lastname@example.org