Organ-specific transcriptome analysis reveals differential gene expression in different castes under natural conditions in Apis cerana

Honeybees are one of the most environmentally important insects, as their pollination of various plant species contributes to the balance among different ecosystems. It has been studied extensively for their unique attribute of forming a caste society. Unlike other insects, honeybees communicate socially by secreting pheromones or by exhibiting specific patterns of motion. In the honeybee industry, the Asian honeybees (Apis cerana) and the Western honeybees (Apis mellifera) are dominant species. However, molecular research on the transcriptomes of A. cerana has not been studied as extensively as those of A. mellifera. Therefore, in this study, caste-specific transcriptional differences were analyzed, which provides a comprehensive analysis of A. cerana. In our dataset, we analyzed gene expression profiles using organs from worker, drone, and queen bees. This gene-expression profile helped us obtain more detailed information related to organ-specific genes, immune response, detoxification mechanisms, venom-specific genes, and ovary development. From our result, we found 4096 transcripts representing different gene-expression pattern in each organ. Our results suggest that caste-specific transcripts of each organ were expressed differently even under natural conditions. These transcriptome-wide analyses provide new insights into A. cerana and that promote honeybee research and conservation.

A comparison of the genomes of A. cerana with two other species of insects, i.e., Drosophila melanogaster and A. mellifera in the NCBI v2.0 database (ACSNU-2.0), provided the coding domain sequence and the read count for gene expression of A. cerana using the Kallisto software. The result showed a total of 10,651 genes. Highthroughput sequencing analysis indicated that worker, drone, and queen bees had differentially expressed genes. After removing the hypothetical protein, ribosomal protein family, and genes having less than 10 transcripts per million (TPM) reads, a total of 4,096 differentially expressed genes were found among all three castes. Raw and processed data are publicly available in the NCBI/GEO database: (http:// www. ncbi. nlm. nih. gov/ geo/) under accession number GSE164333.
Common transcript patterns between worker-drone, worker-queen, and drone-queen. One of the major goals of our study was to obtain high-quality transcriptome data that could be used to predict global changes in the gene expressions of honeybees under natural conditions. Transcriptome profiling showed 269, 196, and 324 genes differentially expressed in the brain, ventriculus, and rectum, respectively, of workers and drones. Venn diagram data showed a total of 269 brain-associated genes in WB and DB, among which 255 genes were only expressed in WB, 11 genes were only expressed in DB, and only three genes were expressed commonly in both WB and DB (Fig. 2a,b). Among the 269 genes in brain, 14 genes were correlated with caste differentiation such as hormone, visual sense, and neuronal signal peptides. Among these genes, ten genes, including odorant receptor (OR) genes were up-regulated in WB. The differential expression of transcripts in the brain indicates that they affect honeybee behaviors 29,43,44 . These results also suggest that workers can sense a wider range of odors types.
Differences in the diets of each caste reportedly affect gut enzymes 45 . We found 39, 78, 238, and 32 genes specifically expressed in WV, DV, WR, and DR, respectively. Among these genes, 30 genes were exclusively expressed in ventriculus, whereas, only five genes were expressed in the rectums of both worker and drone (Fig. 2c). In ventriculus, more transcripts were overexpressed in DV compared to WV (Fig. 2d). In the rectum, however, more transcripts were overexpressed in WR compared to DR (Fig. 2d). The ventriculus and rectum, which are typical digestive organs of honeybees, showed substantial differences in the gene expression profile transcripts, despite being anatomically connected.
The reproductive organs of drones (haploid), workers (diploid), and queen (diploid) were compared. Three commonly expressed genes were found in both DT and DMG, and 18 commonly expressed genes in both WO and QO. The numbers of specific transcripts expressed in DT, DMG, WO, and QO were found in 49, 21, 27, and 418 genes, respectively (Fig. 2e). Transcript profiling between DMG and DT revealed that more genes were www.nature.com/scientificreports/ overexpressed in DT compared to DMG (Fig. 2f). In the comparison between QO and WO as the representative female reproductive organ, most transcripts were overexpressed in QO than in WO (Fig. 2f).
Differentially expressed genes in the brain, gut, and reproductive organs. To understand the honeybee's different caste-specific behaviors several categories of genes related to social communication were investigated. The expression patterns of genes were compared among the worker, drone, and queen. Genes from the brain were classified genes into four categories; OR, neuronal genes, photoreceptors, and hormones ( Fig. 3a and Supplementary Table S2). A comparison of the OR genes between WB and DB showed remarkably high expression profiles in WB. However, it remains to be determined how each OR affects the honeybee's behavior.
In the other two groups (neuronal genes group and photoreceptor group), the genes related with each group were expressed more in DB than in WB. Of the hormone group, prohormone-1 (ACSNU02044T0) was the only hormone gene that was overexpressed in WB. We classified gut genes into three categories digestive enzyme genes, carbohydrate biosynthesis genes, and lipid transport genes. In a previous study of the honey bee protein atlas 19 , major digestive enzymes were highly expressed in workers among castes. Comparing these to our data, we found that transcripts of the digestive enzyme were expressed mostly in the rectum of workers (Fig. 3b). The differential profiles of digestive enzyme genes between the worker and drone suggest that there are differences between the diet composition of the two castes.
To understand the difference between haploid and diploid, we collected and examined the gene expressions of four different reproductive organs of worker, drone and queen (Fig. 3c). Vitellogenin is a precursor protein of egg yolk that is used as a biomarker in female [46][47][48] . In reproductive organs, we found that vitellogenin was expressed in all four reproductive organs. It was most highly expressed in QO, followed by WO, DMG, and DT (Fig. 3c). In a previous data, vitellogenin was first detected in the queen at the mid-late pupal stage, in the worker at the late pupal stage, and in the drone at the adult emergence stage 49 . We classified the genes of reproductive organs into three categories of major royal jelly protein (MRJP) family, hormone-related genes, and growth factors ( Fig. 3c and Supplementary Table S2), which are key regulatory factors during early development. When comparing QO with WO, DMG, and DT, transcripts were mostly expressed in QO (Fig. 3c). To validate the result of RNA-seq data, we randomly selected six genes from the brain and reproductive organs and determined their transcription   www.nature.com/scientificreports/ levels using qRT-PCR. The expression patterns of the six genes based on qRT-PCR showed the same patterns (Fig. 3d), which validated the reliability and reproducibility of the RNA-seq data.
Distribution of the immune system across organ types under natural conditions. To understand the innate immune system of A. cerana, we analyzed the differences in the expression patterns of immunerelated genes between castes and displayed them in a heat map ( Fig. 4a and Supplementary Table S3). We compared 61 immune-related gene expressions in several specific castes. We found that six immune-response peptides (apidaecin, apidermin, hymenoptaecin, hexamerin, defensin, and T-cell immunomodulatory protein) were highly expressed in all the examined organs. It is well-known that apidaecin, apidermin, hymenoptaecin, hexamerin, and defensin are typical AMPs that correspond to the Toll pathway or Imd pathway. Genes related to the Wnt-signal pathway and RNAi pathway showed very low expression patterns of immune-response genes in the reproductive organs of drones and in the guts of workers and drones (Fig. 4a). Among the entire gene expressions of different organs from each caste, we investigated the expression ratio of immune-response genes relative to whole transcripts. In comparing the gene expressions, which are displayed in Fig. 4a, the rate of the expression of immune-related genes was different in each organ, which implies that these genes were mainly abundant in the intestinal organs (Fig. 4b). In comparing the transcriptome profiles of whole organs, queen exhibited weakest innate immune-related genes expressed in natural conditions (Fig. 4c). www.nature.com/scientificreports/ Expression of detoxification machinery. We investigated the expression profile of detoxification machinery in unstressed honeybees. In the wild, forager bees collect pollen and transport them to the colony and nurse bees serve nectar as their nutrient resource to feed the brood. Therefore, it is important to detoxify toxic substances such as pesticides that can contaminate the dietary source of the colony. As the detoxification mechanism is assumed to differ among worker, drone, and queen, we compared the detoxification-related expression patterns of five gene families between different castes under natural conditions. First, there were 12 genes in the glutathione-S-transferase family. In order of high expression ratios were 24% in DT, 19% in WB, and 13% in WR (Fig. 5a). Second, there were 39 genes corresponding to oxidation-related biomolecules (peroxiredoxin, superoxide, superoxide and thioredoxin transcripts). We observed the highest expression ratio in WB was 21% (Fig. 5b). Third, there were two genes in catalase. For which the expression ratios were 23% and 22% in WR and WV, respectively (Fig. 5c). Fourth, there were seven genes in the multidrug resistance-associated protein family. In order of high expression, the ratios were 28% and 24% in WR and QO, respectively (Fig. 5d). Fifth, there were 25 genes in the cytochrome P450 transcripts family. The high expression ratio was 37% in WR, then 22% in WV (Fig. 5e). Genes related to detoxification generally showed high expression rates in WV and WR, which may be closely related to the behavior of the worker acting as nurse, server, and pollinator. Compared to the entire transcripts of each organ, genes encoding peroxiredoxin, superoxide, and thioredoxin were expressed ubiquitously with no organ or caste-specific bias (Fig. 5f), which is a similar pattern to the previous proteome report on A. mellifera 19 .

Differences in the haplodiploid sex-determination system between female ovaries.
To understand haplodiploid development more comprehensively, we compared transcriptome data between fertile QO and unfertile WO. We found a total of 463 genes and displayed them on a heat map (Fig. 6a). These genes were grouped into three clusters. We classified the genes up-regulated in WO as Cluster 1, genes with no difference in expression between QO and WO as Cluster 2, and genes up-regulated in QO as Cluster 3 (Supplementary  Table S4). To obtain information about the biological processes, cellular components, and molecular function of each gene cluster, we analyzed the gene ontology (GO) in three clusters using BLAST2GO program. GO analysis indicated that genes in Cluster 1 were related to the development process, genes in Cluster 2 were related to the development process and immune system process, and genes in Cluster 3 were related to the development process and reproduction (Fig. 6b). Our results showed that both QO and WO were related to the development process, but only QO was related to reproduction process. In molecular biology, transcription factors and elongation factors are key regulators in the transcription of genetic processes 50 . Serine/threonine protein kinase plays a significant role in post-translational modification 51 . www.nature.com/scientificreports/ The proteins that are encoded by these genes are essential for the differentiation and proliferation of germ cells as ovary development factors. In the comparison between QO and WO, QO exhibited extremely large proportions of these factors (Fig. 6c). The transcription levels of nine selected genes encoding for major factors of ovary development can be seen in Fig. 6d. Two of the major factors were up-regulated in WO but others were up-regulated in QO (Fig. 6d). To determine the fold changes in the expression levels of two selected genes of maelstrom and piwi between QO and WO, the results from qRT-PCR analysis matched well with the RNA-seq results (Fig. 6e), which validated the reliability and reproducibility of RNA-seq data.

Expression patterns of venom encoding genes in organs.
There were 105 genes related to venom processes. Initially, high expression patterns were expected for the 105 venom-related genes in WVG, but this was not the case. The expressions of the 105 venom encoding genes in the reproductive organs of the drone and queen were compared with those of the WVG. As expected, genes related to venom protease and phospholipase were generally highly expressed in WVG. However, some venom genes, including phospholipase, prepromelittin gene (which is required to produce melittin), and the venom allergen 5-like (ACSNU08074T0), which is a toxic peptide, showed higher expressions in the drone and queen reproductive organs (Fig. 7a). The mRNA-seq results using qRT-PCR revealed the same patterns for four randomly selected genes (Fig. 7b). These results validate the reliability and reproducibility of RNA-seq data.

Discussion
Here, we presented high-quality transcriptome profiles of organs from different caste members of the A. cerana.
In the data set, our results showed differences in transcripts between female and male honeybees in natural conditions. Furthermore, our transcriptome analysis of A. cerana, provides an important resource for future molecular studies of Asian honeybees, particularly to elucidate the mechanisms of the innate immune system, social behavior, genetic difference between females (worker and queen) and male (drone), and the potential function of haplodiploid development.
AMPs are effective against pathogens that are frequently encountered in honeybees hive 33 . Our data demonstrate that some AMPs were highly expressed in DR, and DMG. The results indicate that AMPs and detoxification-related genes were expressed differently among castes. The expression of AMPs-related genes was mainly expressed, while. detoxification genes were low expression in drones. However, both of the AMPs-related genes and detoxification genes were highly expressed in workers. Therefore, we presume that workers possess immunity against bacteria that are introduced from the outside while foraging outside as a pollinator or as nursing brood to feed pollen. Previous studies showed the workers infected by bacteria had a strong immune response 38 .
It has been reported that the nutrition and dietary habits of bees affect their physiology 52 . Therefore, in the present study, we compared gene expression patterns between worker guts and drone guts. Detoxification mechanisms in honeybee guts, which help bees develop tolerate to a variety of potentially toxic secondary metabolites and pesticides that they encounter in floral nectars and pollen, remain largely unknown 53 . The latest study on the protein atlas of honey bee organs 19 represented the robust detoxification machinery of worker organs. However, www.nature.com/scientificreports/ our expression data shows that the detoxification-related transcriptome expressed more in the drone, than in the worker. In this study, we provide information on these observed transcriptional changes, which are not reflected in data that was obtained on the protein level 19 . The honeybee has a unique sex-determination system that is also known as haplodiploidy. Unfertilized eggs are haploid, whereas fertilized eggs are diploid 54 . Haploid eggs develop into males, whereas diploid eggs develop into females. Vitellogenin is a precursor protein of egg yolk that is used as a biomarker in female [46][47][48] and is also expressed in the fat bodies of workers and queens 55 . In a previous study, the vitellogenin protein played an important role in ovary development, immunity, stress response, and sex-determination [46][47][48]56 . In our results, vitellogenin was expressed in both QO and WO, but transcripts related with ovary development were highly expressed in QO.
Interestingly, the maelstrom transcript and piwi protein transcript were found in both QO and WO. Mael and Piwi proteins are required for transcriptional silencing, which is induced by the piRNA pathway. Our results showed that different expressions of mael and piwi were observed in QO and WO. The Mael protein repressed canonical RNA polymerase II transcription and inhibits germline transposon transcription 57,58 . Based on previous papers and our recent findings, we propose that mael and piwi play an important role in the development of QO and WO.
Among all organs, venom-producing organs give honeybees the ability to protect their colonies from enemies by using venom as a weapon 59 . The WVG produces several components of honeybee venom that cause allergic responses in humans 60 and other vertebrates. Honeybee venom contains the major allergens of Api m 1 (phospholipase A2), Api m 2 (hyaluronidase), Api m 3 (acid phosphatase), Api m 4 (melittin), and Api m 7 (CUB serine protease). Melittin, the main constituent of honeybee venom, is one of the major proteins that composes venom toxicity 61 and is derived from promelittin 62,63 . However, it was not unknown how prepromelittin is converted into promelittin. In our results, prepromelittin was proven to be transcribed in the reproductive organs of drones and queen. Therefore, it is speculated that cleavage occurs with melittin after the gene is transferred and that biosynthesis of melittin from prepromelittin might occur after the prepromelittin gene is transferred from the parent to the worker bee. Our data do not clearly demonstrate this hypothesis, but we provide a clue that indicates the synthesis process of the melittin, which guides future research.
In this study, we investigated the transcriptome profiles in several organs from different castes to determine unique traits. The results revealed that immune-related genes and detoxification mechanisms were generally highly expressed in workers in natural conditions. Thus, we provide a foundation for understanding the physiology, morphology, and caste formation of A. cerana at the molecular level. We hope that this study will stimulate future research on the Asian honeybee, which is less understood compared to Western honeybees.

Methods
Sample collection. Samples from each developmental stage of A. cerana were collected from colonies kept in Cheonan, South Korea. To obtain for the samples, one colony of a mated egg-laying queen was used. Bees were collected from the colony during swarm season, and placed in a 1.5 ml microfuge tube. Prior to dissection, the bees were exposed to CO 2 in order to aid handling. Each organ was dissected using forceps and microdissection scissors. All dissected organs were rinsed with PBS solution 64 . Each dissected organ sample was flash-frozen in liquid nitrogen and then stored at − 80 °C. Organs were obtained from three castes of the A. cerana i.e., worker,  www.nature.com/scientificreports/ drone, queen. Worker organs included worker brain (WB), worker ventriculus (WV), worker rectum (WR), and worker venom gland (WVG). Drone organs include drone brain (DB), drone ventriculus (DV), drone rectum (DR), drone testis (DT), and drone mucus gland (DMG). Queen organ includes queen ovary (QO).
Preparation of the mRNA-sequencing library. Total RNA was extracted from the organ samples using TRIzol reagent (Invitrogen). An equal amount of 5 μg total RNA from each sample was used to construct the mRNA library according to manufacturer instructions (Lexogen, Austria). The libraries were sequenced at a high-throughput sequencing facility (Macrogen, South Korea) on an Illumina HiSeq 2500 platform. The sequencing data were obtained with two biological replicates of each organ.
Computational analysis. The raw Illumina sequence reads were filtered to remove low quality sequences by using fastp 65 prior to bioinformatic analysis. The filtered Illumina short reads were then mapped to the A. cerana reference genome (ACSNU-2.0) 66 using Bowtie2 software 67 . The mapping results of each sample to the reference mRNA sequences obtained by Bowtie2 were then quantified by Kallisto software 68 . The gene expression profiles of each sample obtained by Kallisto were calculated as TPM value to compare the gene expression levels between samples.
Validation of differentially expressed genes by quantitative real-time PCR. The relative mRNA levels were quantified according to the manufacturer's instruction using an AccuPower 2X GreenStar qPCR Master Mix (Bioneer, South Korea). Results were normalized to actin mRNA. Gene expression levels were calculated using the comparative Ct method. All experiments were carried out at least three times for each of three biological replicates. The gene-specific primers used for qRT-PCR are listed in Supplementary Table S5.

Data availability
The datasets generated for this study can be found in the NCBI SRA repository, https:// www. ncbi. nlm. nih. gov/ sra/, with the GEO Accession No.: GSE164333. www.nature.com/scientificreports/