Introduction
The western honey bees (HBs) (Apis mellifera) are social insects and play an important role in honey production, biodiversity maintenance and plant pollination (Potts et al., Reference Potts, Biesmeijer, Kremen, Neumann, Schweiger and Kunin2010). The neurobiological processing signals for crucial nursing behavior are controlled by central, higher order processes in the brain of HBs (Menzel et al., Reference Menzel, Leboulle and Eisenhardt2006). Expression variation in the brain genes are detected for some of the social behavior of that insect, while a large number of HB physiological changes are related to their behavioral plasticity.
Hypopharyngeal and mandibular glands are the two important glands in the head of HBs that their secretion are age dependent and are responsible for secreting the royal jelly (RJ) proteins, diet composition and caste differentiation (Zhang et al., Reference Zhang, Hu, Han, Wei, Meng, Wu, Fang, Feng, Ma and Rueppell2020), response to stresses like starvation and high temperature (Ueno et al., Reference Ueno, Takeuchi, Kawasaki and Kubo2015). RJ is the essential food for broods and queens of HBs (Fujita et al., Reference Fujita, Kozuka-Hata, Ao-Kondo, Kunieda, Oyama and Kubo2013) and contains protein, monosaccharides, fatty acids, vitamins, minerals and water (Viuda-Martos et al., Reference Viuda-Martos, Ruiz-Navajas, Fernández-López and Pérez-Álvarez2008).
The molecular processes involved in age-related social behavior pertain to changes in the brain of adult HBs. For this reason, A. mellifera is the most important species of social insects for studying patterns of gene expression. By sequencing the genome of HB (The Honeybee Genome Sequencing Consortium, 2006) and improving its annotation (Elsik et al., Reference Elsik, Worley, Bennett, Beye, Camara, Childers, de Graaf, Debyser, Deng, Devreese, Elhaik, Evans, Foster, Graur, Guigo, Hoff, Holder, Hudson, Hunt, Jiang, Joshi, Khetani, Kosarev, Kovar, Ma, Maleszka, Moritz, Munoz-Torres, Murphy, Muzny, Newsham, Reese, Robertson, Robinson, Rueppell, Solovyev, Stanke, Stolle, Tsuruda, Van Vaerenbergh, Waterhouse, Weaver, Whitfield, Wu, Zdobnov, Zhang, Zhu and Gibbs2014), several transcriptomic studies carried out to elucidate various aspects of worker HB behavior during development (Hu et al., Reference Hu, Ke, Wang and Zeng2018; Kang et al., Reference Kang, Kim and Lim2021). Studies showed that energy metabolism, catalytic activity, amino acid metabolism, protein synthesis and transport, ribosome pathway and protein processing in the endoplasmic reticulum, TOR signaling pathway, PI3 K-Akt signaling pathway, major royal jelly proteins (MRJP) genes, citrate cycle, fatty acid elongation, protein export, N-Glycan biosynthesis, carbon metabolism, folate biosynthesis, cysteine and methionine metabolism, aminoacyl-tRNA biosynthesis, oxidoreductase activity, ribosome biogenesis, transmembrane transport, energy metabolism, amino acid metabolism, fatty acid, DNA replication. amino acid activation, translation factor activity, and protein modification process differentially contribute in hypopharyngeal glands (HPGs) development (Liu et al., Reference Liu, Ji, Yin, Shen, Shen and Chen2013; Nie et al., Reference Nie, Gao and Zhu2021). Also, genes necessary for synaptic/neurotransmission, receptor signaling pathways, nervous system development, protein folding G-protein coupled receptor signaling, protein kinase activity, insulin receptor signaling, and response to heat that workers needs for making spatiotemporal memories of foraging are differentially expressed in the head of nurse and forager HBs (Khamis et al., Reference Khamis, Hamilton and Medvedeva2015).
Despite the major research focuses in the field of genomic or transcriptomic studies of HB, a few studies have been conducted on encoding genes and pathways influenced by head development in worker HBs (Drapeau et al., Reference Drapeau, Albert, Kucharski, Prusko and Maleszka2006; Buttstedt et al., Reference Buttstedt, Moritz and Erler2014; Nie et al., Reference Nie, Liu, Pan, Li, Li, Zhang, Chen, Miao, Zheng and Su2017; Dobritzsch et al., Reference Dobritzsch, Aumer, Fuszard, Erler and Buttstedt2019). The aim of this research was to study the transcriptome of head capsule of HB workers through comparing the gene expressions at 1 and 9 days after eclosion by focusing on RJ production and nursing behavior.
Material and methods
Material preparation
Apis mellifera meda is selected for the current study because it is a native subspecies that dispersed throughout Iran and also has adapted to undesirable environmental conditions of Iran (Ebadi and Ahmadi, Reference Ebadi and Ahmadi2016; Modaber et al., Reference Modaber, Nazemi-Rafie and Rajabi-Maham2019). Colony selection was carried out based on Modaber et al. (Reference Modaber, Nazemi-Rafie and Rajabi-Maham2019). The colonies of A. m. meda were supplied from the apiary of University of Kurdistan, Iran. Ten mated HB queens of the same age were inserted randomly into ten hives as one queen existed in each hive. The combs containing the last pupal stage were transferred into screened cages. They placed into an incubator with 75–80% relative humidity and a temperature at 32°C until metamorphosis into the adult workers. Newly emerged HBs (<1 h) were marked using white and red colors and transferred into the parent hives (fig. 1). Then, the red and white marked workers were gathered on nine days and one day after eclosion (the emergence of a worker from the pupa), respectively. 100 1_DO HBs and 100 9_DO HBs selected from each hive. The workers were stored at − 80°C prior to dissection.
Marked worker HBs at 1_DO and 9_DO. (a) 9_DO worker HBs, (b) 1_DO workers HBs.
Sample preparation, RNA extraction and sequencing
The heads of selected workers of each hive were detached and pooled for RNA extraction in liquid nitrogen. RNA extractions were performed by RNeasy Mini Kit-QIAGEN following the manufacturer's protocol. Afterwards the concentration, A260/230 and A260/280 of extracted RNA were assessed by a nanodrop. In addition, gel electrophoresis was used for checking the RNA integrities. Total RNA sequencing of RNA samples were done with an Illumina HiSeq system (Genewiz company, USA).
Analysis of gene expression
Quality control and trimming of sequenced reads were done by FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and Trimmomatic (Bolger et al., Reference Bolger, Lohse and Usadel2014), respectively. Then reads mapped to the Amel_3.5 assembly by using the Tophat2 (Kim et al., Reference Kim, Pertea, Trapnell, Pimentel, Kelley and Salzberg2013) and the gene expressions quantified by HTSeq (Anders et al., Reference Anders, Pyl and Huber2015). DESeq package were used for differential gene expression and functional analysis were carried out with DAVID and KEGG databases.
PPI network and hub gene identification
Interacting differential expressed genes (DEGs) were retrieved from STRING network (https://string-db.org/; v11.0) (Szklarczyk et al., Reference Szklarczyk, Gable, Lyon, Junge, Wyder, Huerta-Cepas, Simonovic, Doncheva, Morris and Bork2019) for construction the protein-protein interaction networks (PPI). The PPI visualized by Cytoscape (v3.8.2) (Shannon et al., Reference Shannon, Markiel, Ozier, Baliga, Wang, Ramage, Amin, Schwikowski and Ideker2003). CytoHubba, based on the maximal Clique Centrality (MCC) algorithm plugin, was used for identification of the hub genes in the co-expression network of the current study. The MCC, MNC, Degree, EPC and EcCentricity algorithms of this plugin were employed for detecting the significant genes. Indeed, the PPI networks of selected important genes included in fatty acid and protein related pathways, heat shock proteins and histone-lysine N-methyltransferase encoding genes were identified as smaller networks. The Degree, one of the centrality parameter of CentiScape plugin of Cytoscape, employed for finding the significant gene of small PPI networks.
Results
RNA-sequencing
Sequencing the RNA samples generated an average of 55 million reads of paired end sequences with 150 bp length. 41458895 and 41006226 reads of 1_DO and 9_DO workers, respectively, aligned to A. mellifera reference genome.
Analysis of gene expression
The results showed that 1662 gene transcripts have differentially significant expressions between 9_DO and 1_DO HBs. Of these, 1125 had higher expressions in 9_DO HBs and 537 genes up regulated in 1_DO (fig. 2 and tables 1 and 2). As shown in tables 1 and 2, some biological pathways with significant expressions were common between the two studied groups and differed only in term of active genes in their pathways. But in some cases, certain pathways differentially expressed in each group. In the following, important pathways with differential expressions are considering in details.
Volcano plot of head gene expression pattern between 9 and 1 day old HBs. Gray dots indicate the non-DEGs. The green and red dots represent the down regulated and up regulated genes (P < 0.05) in the head of 9_DO and 1_DO HBs, respectively.
List of up-regulated KEGG pathways and their related IDs in 9_DO HBs compared with 1_DO HBs
*The underlined pathways are specifically up-regulated in 9_DO HBs.
List of down-regulated KEGG pathways and their related IDs in 9_DO HBs compared with 1_DO HBs.
*The underlined pathways are specifically down-regulated in 9_DO HBs.
Major royal jelly proteins (MRJP)
Expression of genes encoding MRJP 1, 2, 3, 5, 6 and 9 was significantly higher in 1_DO HBs (fig. 3a). As shown in tables 1 and 2 and also table S1, the biosynthesis and metabolism of various amino acids and protein processing in endoplasmic reticulum pathway are statistically more active in 9_DO HBs compared with 1_DO HBs.
Comparisons of differential gene expression as heat map. Heat map represented the significant gene expression changes of (a) RJ protein, (b) Fatty acid related pathways (c) Hippo signaling pathway.
Fatty acid metabolism
Most of the differentially expressed genes relevance to fatty acid pathways was up-regulated in 9_DO HBs. They are glycosphingolipid biosynthesis, fatty acid metabolism, fatty acid biosynthesis, fatty acid elongation, fatty acid degradation, glycerolipid metabolism, glycerophospholipid metabolism, ether lipid metabolism, sphingolipid metabolism, glycosphingolipid biosynthesis and biosynthesis of unsaturated fatty acids (table 1). As shown in fig. 3b, in 1_DO HBs, only fatty acid degradation pathway had higher expression as a result of higher activity of alcohol dehydrogenase class-3 and glutaryl-CoA dehydrogenase genes, while the number and expression levels of genes contributed in fatty acid biosynthesis, elongation and metabolism are extremely higher in 9_DO HBs (tables S1 and S2).
Starch and sucrose metabolism
Transcripts associated with Starch and sucrose metabolism with high expression in 9_DO HBs were alpha-amylase, Hbg2, Hbg3 and trehalase. However, the glycogen synthase, glycogen debranching enzyme, hexokinase type 2, hexokinase-1-like and trehalase-like genes were down-regulated in this group (table S2).
Hippo signaling pathway
This is a conserved pathway that controls HB organ size by regulating cell proliferation and cell death processes. It is interesting that the expression of genes encoding casein kinase I-like, partitioning defective 3 homolog, protein kibra, and scaffold protein salvador-like in 9_DO HBs decreased compared to 1_DO HBs. While, the expression of 14-3-3 protein zeta, LIM domain-containing protein jub, ras association domain-containing protein 2, serine/threonine-protein phosphatase 2A catalytic subunit alpha isoform, and transcriptional enhancer factor TEF-1 genes were higher in 9_DO HBs.
Transforming growth factor beta 1 (TGF-β) induced pathways
This contributed in wound-healing and plugging responses in HBs (Richards et al., Reference Richards, Jones and Bowman2011), is down regulated in 9_DO HBs by reduction in the expression of cullin-1 and dorsal-ventral patterning protein Sog genes (table 1 and tables S1 and S2). On the other hand, Jak-STAT signaling, Wnt signaling, mTOR signaling, and MAPK signaling pathways are up regulated by higher expressions of G1/S-specific cyclin-D2, signal transducing adapter molecule 1, G1/S-specific cyclin-D2, 3-phosphoinositide-dependent protein kinase 1, RAC serine/threonine-protein kinase, mitogen-activated protein kinase 1, serine/threonine-protein kinase mTOR, serine/threonine-protein kinase unc-51, and target of rapamycin complex subunit lst8 genes. In addition, the expression of toll pathways related genes as toll, toll-interacting protein, and tolloid-like protein 1 increased in 9_DO bees (tables S1 and S2). It is necessary to mention that high number of genes that contribute to phagosomes and phagocytosis like Mig-2-like GTPase Mtl, V-type proton ATPase subunit B, V-type proton ATPase subunit d, lysosome-associated membrane glycoprotein 1, protein transport protein Sec61 subunit alpha, protein transport protein Sec61 subunit beta, tubulin beta chain-like, vacuolar H + ATP synthase 16 kDa proteolipid subunit, ras-related protein Rab-7a, ras-related protein Rac1, and etc. up-regulated in 9_DO HBs (tables S1 and S2).
Cytochrome P450 (CYP450) genes
Members of this group have a diverse expression during HB development. The results exhibit greater expression of cytochrome P450 305a1, cytochrome P450 4C1, cytochrome P450 6k1 and cytochrome P450 4c3 in 9_DO HBs while, probable cytochrome P450 6a14, cytochrome P450 18a1, cytochrome P450 6a17 and cytochrome P450 304a1 down-regulated in the older HBs (tables S1 and S2).
39S ribosomal genes
Most of the these genes (L12, L13, L30, L37, L44, L46, L53, L9, S24 and L24), ribosomal RNA processing protein 1 homolog and ribosomal protein S6 kinase beta-1-like had higher expressions in 9_DO HBs. Except for 40S ribosomal protein S24, the gene expression of the other 40S ribosomal proteins (S11, S12, S15Aa, S16, S2, S21, S3a, S7) and 60S ribosomal proteins (P0, P2, 13a, L15, L17, L22, L29, L37, L44, L7, L8 and L9) were higher in 1_DO HBs (table S2).
Vitamin B6 metabolism pathway
In the present study, this pathway had a higher expression in 9_DO HBs through the differentially transcription process of pyridoxal kinase and pyridoxine-5′-phosphate oxidase-like genes (table S2).
RNA degradation pathways
Most genes encoding proteins active in this pathways, including 60 kDa heat shock protein, activator of 90 kDa heat shock protein ATPase homolog 1, heat shock protein cognate 5, dnaJ homolog subfamily C member9-like /3/ 16/ 11 and dnaJ homolog subfamily B member 11/ 12 had higher expressions in 9_DO HBs. By contrast, the expression of Hsc70-3, heat shock protein beta-1, heat shock protein 83, heat shock protein Hsp70Ab-like, Hsp90, dnaJ homolog subfamily C member 22 and dnaJ protein homolog 1 encoding genes were down regulated in 9_DO HBs compared to the 1_DO (tables S1 and S2).
Other pathways
The results of this study showed that the activity of Phenylalanine, tyrosine and tryptophan biosynthesis pathway reduced in 9_DO HBscompared with 1_DO (table 2). That occurred by the decrease in expression of the tyrosine aminotransferase gene which is involved in the production of 4-hydroxy-phenylpyruvate, phenylpyruvate, tyrosine, phenylalanine and ultimately glucosinolate biosynthesis.
It is obvious that neuro-transcriptome of HB strongly connected with the bee social characteristics. Neuroactive ligand-receptor interaction is another pathway that down regulated in 9_DO HBs (table 2). The gene encoding partitioning defective 3 homolog protein had a lower expression in older bees (table S1).
In 9_DO HBs, the expression of eight genes encoding 26S proteasome non-ATPase regulatory subunits was significantly higher (table S1). However, none of the genes active in this process was significantly expressed in 1_DO HBs.
Most of the genes with differential expression and their related pathways especially monosaccharide biosynthesis, pentose phosphate pathway, amino sugar and nucleotide sugar metabolism, citrate cycle, nitrogen metabolism, metabolic pathways, homologous recombination, RNA degradation, carbon metabolism, pantothenate and CoA biosynthesis, amino sugar and nucleotide sugar metabolism, and mismatch repair were up-regulated in 9_DO HBs. They participate in energy supplying, transcription correct or hydrocarbon storage (tables 1 and 2, table S1). Because of they are extremely enrolled in every kinds of head cell, we cannot connect them with professional behavior or activity.
PPI networks
We found a network with 1748 and 8844 nodes and edges, respectively, with 5.06e-10 PPI enrichment p-value from searching in STRING database. These include 55 local network cluster, 27 KEGG pathways, and 6 reference publications. After removing a small number of single proteins out of main networks, we found a network with 1346 nodes in Cytoscape software. As shown in fig. 4, most of the DEGs had higher expression in 9_DO HBs in comparison with 1_DO HBs. The genes that have more interactions are generally that have relatively high or low differential expression. Genes with very high/low expressions have little relation with other genes in the protein network. In studying the PPI network of important genes of fatty acid related pathways, heat shock proteins and histones, we found the followings:
PPI networks obtained from DEGs. The ratio of expressed genes (9_DO HBs /1_DO HBs) is increased from blue to red. Each node represents a protein, and each edge represents an interaction.
Fatty acid related pathways
The first differential expressed gene in fatty acid related pathways was trifunctional enzyme subunit beta that had high interactions with other differential genes. This gene interacted with Hmgs, carnitine O-palmitoyltransferase 2, ATP-citrate synthase, short-chain-enoyl-CoA hydratase-like, dolichyl-diphosphooligosaccharide, acyl-CoA dehydrogenase, uncharacterized (LOC551958), acyl-CoA synthetase, Idh, glycogen synthase, trifunctional enzyme subunit beta, alpha-keto acid dehydrogenase, 3-hydroxyacyl-CoA dehydrogenase type-2-like, SNF1A, 3-ketoacyl-CoA thiolase, citrate synthase 2, Etfdh, cytochrome b-c1 complex, dihydrolipoyllysine, peroxisomal targeting signal 1 receptor, lipid-transfer, electron transfer flavoprotein, acyl-CoA dehydrogenase, ATP5G2, trans-2-enoyl-CoA reductase, and succinyl-CoA:3-ketoacid coenzyme A transferase (tables S2 and S3). As shown in fig. 5a, acyl-CoA synthetase, 3-hydroxyacyl-CoA dehydrogenase, 5′-AMP-activated protein kinase and cytochrome b-c1are down regulated in 9_DO, while the other up regulated. α-keto acid dehydrogenase complex is responsible to α -keto acids oxidative decarboxylation for entering into TCA cycle or lipid biosynthesis (Nobukuni et al., Reference Nobukuni, Mitsubuchi, Akaboshi, Indo, Endo, Yoshioka and Matsuda1991). Acyl-CoA synthetase is the activator of fatty acid metabolisms (Watkins et al., Reference Watkins, Maiguel, Jia and Pevsner2007) and AMP-activated protein kinase is a regulator of energy state in cells (Hardie and Ashford, Reference Hardie and Ashford2014; Nguyen, Reference Nguyen2017) and by increasing that level, it reduces the ATP producing process through turning off the pathways such as the fatty acids and glycogen biosynthesis (Hardie, Reference Hardie2011). The evidence reviewed here seems to suggest that this co-expression network is one of the good samples of controlling interaction in the cells, as by growing the HB cells from 9_DO HBs to 1_DO HBs, the energy producing activity of cells increased and cells uses control processes to adjust the energy balance.
PPI networks of important differential genes involved in fatty acid related pathways and histone encoding genes. PPI networks of (a) Trifunctional enzyme subunit beta, mitochondrial (LOC412345); (b) GB41390-PA (Glycerol-3-phosphate dehydrogenase); (c) Putative glucosylceramidase 4 (LOC409709); (d) AChE-2 (Acetylcholinesterase 2); (e) PR-Set7 (histone-lysine N-methyltransferase). The ratio of expressed genes (9_DO HBs /1_DO HBs) is increased from blue to red. Each node represents a protein, and each edge represents an interaction.
Glycerol-3-phosphate dehydrogenase which contributes in glycerophospholipid metabolism pathway is one the gene with significant expression in 9_DO HBs. It has interactions with glycerol-3-phosphate acyltransferase 1, lysophospholipid acyltransferase 2, Pgk, uncharacterized LOC413484, 1-acyl-sn-glycerol-3-phosphate acyltransferase alpha-like, adenylyltransferase and sulfurtransferase, glycerol-3-phosphate acyltransferase 4, Dhap-at, lysocardiolipin acyltransferase 1-like, 1-acyl-sn-glycerol-3-phosphate acyltransferase alpha-like genes (fig. 5b, tables S2 and S3). All the genes except the later have higher expression in 9_DO HBs. All of them are participating in glycerophospholipid, glycerolipid metabolic metabolism process and have an internal pathway interactions excluding LOC413484 which is an exosome complex rna-binding protein.
Putative glucosylceramidase 4 is one of the gene involved in fatty acid pathways in lysosomes, has a higher expression in 9_DO HBs. Results showed that its protein interacts with uncharacterized LOC410793, Kr, sorbitol dehydrogenase-like, cyclin-G-associated kinase, E3 ubiquitin-protein ligase parkin in HBs (fig. 5c). The last 3 proteins are involved in Protein processing in endoplasmic reticulum, Polyketide synthase, enoylreductase and ATP binding, respectively. Glucosylceramidase 4 has inverse interaction with Kr and cyclin-G-associated kinase as by increasing its expression, the expression of these two genes increased, while it has direct interactions with other mentioned genes (tables S2 and S3). Polyketide synthase which syntheses the aromatic polyketides with highly diverse polyketide chains (Wang et al., Reference Wang, Zhang, Chen, Sun, Yan, Shen and Yuan2020) and also host-associated microbes and biosynthesis of metabolites with antifungal activity (Miller et al., Reference Miller, Smith and Newton2020). This further supports the results of higher expression of immune related genes in 9_DO of the current study.
AChE-2 is the other gene with lower expression in the Glycerophospholipid metabolism of 9_DO HBs. This gene with AChE acitivity interacts with vesicular acetylcholine transporter and CYP18A1 genes (fig. 5d, tables S2 and S3). The CYP18A1 that involved in insect hormone biosynthesis pathway, its expression reduced in 9_DO HBs, while the expression of vesicular acetylcholine transporter increased. In HBs, cytochrome P450 acts in detoxification process of phyto-chemicals in honey and pollen, mycotoxins in the hive environment, pyrethroid pesticides and in-hive acaricides (Mao et al., Reference Mao, Rupasinghe, Johnson, Zangerl, Schuler and Berenbaum2009, Reference Mao, Schuler and Berenbaum2011, Reference Mao, Schuler and Berenbaum2017; Niu et al., Reference Niu, Johnson and Berenbaum2011). Also, the expression of AChE in response to stress situation like brood rearing suppression, neonicotinoids, crowding and heat shock that faced worker bees is increased (Boily et al., Reference Boily, Sarrasin, DeBlois, Aras and Chagnon2013; Kim et al., Reference Kim, Kim, Lee, Han and Lee2019). The duty of vesicular acetylcholine transporter is carrying acetylcholine into synaptic vesicles (Giboureau et al., Reference Giboureau, Aumann K, Beinat and Kassiou2010) and its role is essential for development (de Castro et al., Reference de Castro, De Jaeger, Martins-Silva, Lima, Amaral, Menezes, Lima, Neves, Pires and Gould2009).The study of this co-expression network shows that the low gene expression of AChE-2 and cytochrome P450 and increased expression of basic gene involved in development is due to the absence of unknown stress and continuing the grow and development of bees from 1 to 9 days.
Histone encoding genes
The pr-set7 is the other important gene with higher expression in 9_DO HBs. In the network formed based on this gene, the eukaryotic translation initiation factor 3 subunit K, DNA damage-binding protein 1, cullin-4A, transforming growth factor beta regulator 1, sGC-alpha1, Cdc6, mus209, lysine-specific histone demethylase 1A (Su(var)3-3), CTR9 and calmodulin-lysine N-methyltransferase are detected as interacting genes which last fourth gene have higher expression in 1_DO HBs (fig. 5e, tables S2 and S3). Considering the result of studying this network, higher expression of transforming growth factor beta regulator, translation initiation factor 3, Cdc6 genes and their accelerating role in cell growth, suggest that this network has a role in increase the cell growth. The Degree score of main network showed that mus209, Cdc6 and translation initiation factor 3 can introduce as key genes of this network (table S2).
Heat shock proteins
The Hsc70-5 is the first heat shock protein member that its PPI network studied in this section. Its PPI network and node table are presented in fig. 6a and table S3.The number of interacting genes with lower expressions in this network is relatively high. What attracts attention in this network in addition to three ribosomal proteins RpLP0, RpS2 and RpL8, two heat shock protein Hsp70Ab and 83 is that they are among the interactive groups with low expression. Indeed, translation elongation factor 2, Atp5a1and Pgk are the three hub genes of this network with the highest Degree core of main network (tables S2 and S3).
PPI networks of important differential heat shock encoding genes. PPI networks of (a) Heat shock protein cognate 5 (Hsc70-5); (b) Heat shock protein 75 kDa (Trap1); (c) 60 kDa heat shock protein (LOC409384); (d) Heat shock protein 83 (LOC411700). The ratio of expressed genes (9_DO HBs /1_DO HBs) is increased from blue to red. Each node represents a protein, and each edge represents an interaction.
The second heat shock protein encoding gene with significantly higher expression is Trap1. What is clear about the Trap1 related network is that its network is simpler than other selected heat shock protein (fig. 6b). Therefore, we can study its components with more details. Cdc37, Hsc70-3 and Hsp70Ab are the interactive genes with lower expression (tables S2 and S3). Crc is not expressed in our samples and it recognized via string interaction search. LOC552527, Hsc70-5, Akt1, Atp5b, 60 kDa heat shock protein, mRpL12, LOC412409, Ndufa10, bor, PpD3, wee, E3 ubiquitin-protein ligase parkin, activator of 90 kDa heat shock protein ATPase homolog 1, CHORD, and crc are the Trap1 interacting genes with higher expression. Overall, they involved in Oxidative phosphorylation, Metabolic pathways, RNA degradation, ATP-binding, Nucleotide-binding, Oxidative phosphorylation, Metabolic pathways, Ubiquitin mediated proteolysis, Protein processing in endoplasmic reticulum, ATP binding, proton-transporting, ATP synthase activity, ATPase activator activity, Jak-STAT signaling pathway (tables S2 and S3). Analysis the different parameters of this network, based on the degree score of main network, revealed that the LOC552527, Hsc70-3 and Hsc70-5 are the three candidate hub genes (tables S2 and S3).
60 kDa heat shock protein, mitochondrial-like is one of the important heat shock protein member with higher expression. This protein interacts with 48 proteins that most of them have higher expression and contribute in Ion transport, ATP synthesis, Electron transport, ATP-binding, TCA cycle, Glyoxylate and dicarboxylate metabolism, Metabolic pathways, amino acids metabolisms (fig. 6c). LOC408734, Hsp90, Hsc70-5 are active in protein folding and response to stress in addition to RpL8, mRpL4 which are involved in translation, have reverse interaction with the 60 kDa heat shock protein. dnaJ homolog subfamily C member 22 with protein refolding function (da Silva Menegasso et al., Reference da Silva Menegasso, Pratavieira, da Gama Fischer, Carvalho, Roat, Malaspina and Palma2017), also down-regulated (tables S2 and S3). The three putative hub genes of this network, based the main network scores, are cell division cycle 5-like protein, Atp5a1 and Hsc70-3.
The heat shock protein 83 is one the selected heat shock protein for PPI network studies. The difference of this gene with other member of this group is that its expression in 9_DO HBs is reduced. Only one encoding 40S ribosomal protein S3a gene was identified in this network with higher expression in 1_DO HBs. Hsc70-3, and heat shock protein Hsp70Ab have lower expression (fig. 6d and tables S2 and S3). What is noticeable in the study of this network is that interacting genes have a high diversity in their activities compared to the other heat shock protein networks (fig. 6). Based on the Degree score of main network, cell division cycle 5-like protein, translation elongation factor 2 and Atp5a1 genes can introduce as the three identified hub genes of this network. Overall, it seems that heat shock protein 83 is a part of another general process instead of being the center of a specific process in 9_DO HBs.
Identified hub genes
The results of identifying 20 hub genes using MCC, MNC, Degree, EPC, and EcCentricity algorithms of CytoHubba plugin are given in table 3. There was a great deal of variation in the results of the used algorithms in terms of the identified genes. The common genes identified between the different methods were shown in the van diagram in fig. 7. We cannot identify any common gene in all applied algorithms. It seems this difference is due to the number of DEGs and variation in their expressions. To avoid confusion, the results of MCC algorithm, most widely used method, are discussed in detail only.
The blue, red, light orange, dark orange, green and pink areas represent the hub genes identified by MCC, MNC, Degree, EPC, and EcCentricity methods. Common parts represent the prevalent hub genes identified by those methods and the unique correspond to uncommon identified genes.
20 top identified genes from PPI network by MCC, MNC, Degree, EPC, and EcCentricity algorithms of CytoHubba plugin of Cytoscape software.
* and ^ represent the up-regulated and down-regulated genes in 9_DO HBs compared with 1_DO HBs, respectively.
Gene name guidance:.
TER94, transitional endoplasmic reticulum ATPase; Hsc70-3, heat shock 70 kDa protein cognate 3; 107964351, 60S ribosomal protein L7; TBL3, transducin beta-like protein 3; GB45948-PA, pseudouridylate synthase 7 homolog; Pgk, phosphoglycerate kinase; Smid, nuclear valosin-containing protein-like; GB42997-PA, uncharacterized; RpS16:40S ribosomal protein S16; LOC551689, nucleolar protein6; Sec61alpha, protein transport protein Sec61 subunit alpha; Atp5a1, ATP synthase subunit alpha, mitochondrial; RpL9, 60S ribosomal protein L9; RpS3A, 40S ribosomal protein S3a; CYP305D1, probable cytochrome P450 305a1; Pix, ATP-binding cassette sub-family E member 1; RpS2, 40S ribosomal protein S2; GB53954-PA, nucleolar protein 56; RpI135, DNA-directed RNA polymerase I subunit RPA2; RpLP0, 60S acidic ribosomal protein P0; RpL8, 60S ribosomal protein L8; LOC726350, putative ATP-dependent RNA helicase DHX33; LOC727591, RRP12-like protein; OC726961, testis-expressed sequence 10 protein homolog; LOC725618:myb-binding protein 1A; LOC552664, probable ATP-dependent RNA helicase DDX28; LOC552527, cell division cycle 5-like protein; LOC412817, importin-4-like; LOC552509, MKI67 FHA domain-interacting nucleolar phosphoprotein-like; LOC408687:putative tRNA (cytidine(32)/guanosine(34)-2′-O)-methyltransferase; LOC410414, 116 kDa U5 small nuclear ribonucleoprotein component; LOC409866, eukaryotic initiation factor 4A-III; LOC552439, WD repeat-containing protein 3; LOC551824, putative pre-mRNA-splicing factor ATP-dependent RNA helicase PRP1; LOC551781, pre-rRNA-processing protein TSR1 homolog; LOC551125, 40S ribosomal protein S15Aa; LOC409167, translation elongation fac; LOC551087, periodic tryptophan protein 2 homolog; LOC412931, U3 small nucleolar RNA-associated protein 15 homolog; LOC412221:proliferation-associated protein 2G4; LOC412190, probable ATP-dependent RNA helicase DDX47; LOC411510, DNA-directed RNA polymerases I and III subunit RPAC1; LOC411250, DEAD-box ATP-dependent RNA helicase 20; LOC411091, RNA-binding protein 34-like; LOC411026, WD repeat-containing protein 36; LOC411024, DDB1- and CUL4-associated factor 1; LOC410799, nucleolar GTP-binding protein 1.
Discussion
Considering the expression pattern of genes in HB heads showed that there were significant expression changes in the wide range of genes between two groups. It indicates that the HBs undergo significant expression changes in their transcriptomes from the exit through the larval stage compared to the RJ production stage. Due to the high number of identified genes and for better understanding, genes are discussed by contributed pathways in this section.
Proteins and amino acids related pathways
RJ is a bee-specific protein that, playing a role in the determination of the developmental fate of larva (Buttstedt et al., Reference Buttstedt, Moritz and Erler2014), is used as a marker for checking the authenticity and quality of honey (Bilikova et al., Reference Bilikova, Krakova, Yamaguchi and Yamaguchi2015) and the base material of some drugs. Due to the fact that 82–90% of the total RJ is composed of MRJP proteins (Schmitzova et al., Reference Schmitzova, Klaudiny, Albert, Schröder, Schreckengost, Hanes, Judova and Šimúth1998), therefore, studying their gene expressions has a high value in the field of RJ production. Dobritzsch et al. (Reference Dobritzsch, Aumer, Fuszard, Erler and Buttstedt2019) showed that the expression of RJ proteins is influenced by genes and age. Indeed, Feng et al. (Reference Feng, Fang and Li2009) by comparing Italian bees with RJ bees reported that RJ HBs can secrete RJ on third day whereas the other do that on day 6. In the present study, the expression level of genes encoding this protein in 1_DO HBs was higher than 9_DO HBs and it is in accord with previous researches. Due to the fact that the amount of RJ in 9_DO HBs reached into the maximum, it seems that it is necessary for the cells to increase the expression level of the relevant genes a few days before (1th day) to achieve the highest amount of RJ on the 9th day. According to these data, we can infer that first day of nurse HB's life is a key stage in the expression of RJ related genes. Therefore, in order to achieve the potential of HB genotypes, producers recommended making effort in order to minimize the occurrence of light, temperature and even food stress. In addition, although highest production of RJ occurs in 9th day of HB life, this step has no significant effect on the expression of genes associated with this protein as their expression decreased over the time. Instead, 9_DO HBs focus on keeping activity of pathways associated with the production and export of this protein.
Because genes contributing in biosynthesis and metabolism of amino acids have important role in various cell progresses, the detailed study of these genes will not yield a definite result. Overall, we can conclude that the protein and amino acid related pathways including biosynthesis, metabolism, transport, degradation or regulation pathways are statistically more active in 9_DO HBs, this can be related to increase the bee activity, cell needs or RJ production and secretion.
The most surprising aspect of the data is related to genes contributed in protein processing in endoplasmic reticulum pathway, that 9_DO HBs differentially overcome 1_DO HBs in terms of number and of expression level of these genes. These results are in agreement with findings which showed the expansion of rough endoplasmic reticulum start from the nurse emergence and stop in foragers (Smodiš Škerl and Gregorc, Reference Smodiš Škerl and Gregorc2015).
Fatty acid related pathways
The focus of 9_DO HB is on the fatty acid elongation and degradation process rather than fatty acid biosynthesis. As only fatty acid synthase gene with biosynthesis activity has higher expression in this group but five genes (3-ketoacyl-CoA thiolase, trifunctional enzyme subunits alpha and beta, very-long-chain 3-oxoacyl-CoA reductase and very-long-chain enoyl-CoA reductase) and nine genes (3-ketoacyl-CoA thiolase, acetyl-CoA acetyltransferase, alcohol dehydrogenase class-3, carnitine O-palmitoyltransferase 2, medium-chain specific acyl-CoA dehydrogenase, trifunctional enzyme subunits alpha and beta, very long-chain specific acyl-CoA dehydrogenase, mitochondrial, and retinal dehydrogenase 1) had higher expressions in elongation and degradation process (table S2). In worker HBs, mandibular glands play an important role to feed larvae in active period of the HPGs, because they produce the fatty acids of 10-HDA and 10-HDAA (Winston, Reference Winston1991; Isidorov et al., Reference Isidorov, Bakier and Grzech2012; Li et al., Reference Li, Huang and Xue2013). 10-HAD and 10-HDAA are the most important fatty acids of mandibular glands (60–80% of the total fatty acid composition of RJ) that are added to RJ produced from HPGs (Keeling et al., Reference Keeling, Otis, Hadisoesilo and Slessor2001; Isidorov et al., Reference Isidorov, Bakier and Grzech2012; Li et al., Reference Li, Huang and Xue2013). The secretions of fatty acids of mandibular glands coincide with active period of the HPGs (Keeling et al., Reference Keeling, Otis, Hadisoesilo and Slessor2001). The HPGs are inactive in 1_DO HBs therefore mandibular glands are not active too. In this research, the metabolism of fatty acids differentially increased in 9_DO HBs therefore higher expressions in elongation and degradation process were found in 9_DO HBs compared with 1_DO HBs.
It seems that the biosynthesis and metabolism of fatty acids are not a priority for 1_DO HBs and they implied most of the produced lipids for RJ production or partial developmental progress. Indeed the metabolism of fatty acids differentially increased by HB growth up to day 9th.
Starch and sucrose metabolism
Among the DEGs contribute in the Starch and sucrose metabolism, glycogen debranching enzyme involves in Glycogen degradation that its down-regulation in 9_DO HBs help to improve the glucose level by catalyzing the breakdown of glycogen. The up regulated of alpha-amylase genes is due to role of HPGs in honey production by producing converting enzymes for making sucrose from nectar. This result is in agreement with the finding of Vannette et al. (Reference Vannette, Mohamed and Johnson2015) that reported by increasing the age of HBs, the alpha-amylase transcripts in HPGs transcriptome enriched in foragers.
Immune system related pathways
The 18-wheeler coding gene had a higher expression in 9_DO HBs, which due to its role in melanogaster morphogenesis (Eldon et al., Reference Eldon, Kooyer, D'Evelyn, Duman, Lawinger, Botas and Bellen1994) and bee immune system (Riddell et al., Reference Riddell, Garces, Adams, Barribeau, Twell and Mallon2014; Fine et al., Reference Fine, Cox-Foster and Mullin2017). It seems that with increasing bee life duration, its immune system has been strengthened by increasing the expression of this gene.
The HBs induce their immune responses via cellular response and humoral pathways. Phagocytosis, encapsulation and nodulation process are the cellular responses and humoral immunity consists of antibody-mediated immunity and phenoloxidase activity (Hystad et al., Reference Hystad, Salmela, Amdam and Münch2017). Studies have shown that the number of phagocytes in forger bees significantly increases compared to nurse bees due to the increase in the number of dead cells (Amdam et al., Reference Amdam, Aase, Seehuus, Fondrk, Norberg and Hartfelder2005). The results of present study which showed the expression of phagosome related genes differentially increased from 1 to 9_DO HBs further support the idea of previous studies.
Toll pathway includes Toll and the Toll-like receptors are belonged to transmembrane signal transducing proteins and act as immunity and development factors (Evans et al., Reference Evans, Aronstein, Chen, Hetru, Imler, Jiang, Kanost, Thompson, Zou and Hultmark2006). Erban et al. (Reference Erban, Sopko, Kadlikova, Talacko and Harant2019) by comparing the Varroa and non-parasitized HB proteomes reported that Varroa activates TGF-β-induced, Jak-STAT signaling, MAPK cascade, and mTOR signaling pathways. There are two possible explanations for the statistically significant differences in mentioned pathways between 1 and 9_DO HBs as infection of HBs with pathogens or bee planning to strengthen the immune system to prepare to leave the hives. As respect that no pathogenic factor or special stress observed in the hives of this experiment, it can be concluded that the natural growth cycle of bees is planned to produce gradually proteins involved in the bee immune system or it can be due to the hazardous microorganisms and chemical compounds that entered by forager bees into the hives. These results are in accord with Vannette et al. (Reference Vannette, Mohamed and Johnson2015) indicating that the genes related to immune signaling pathways had higher expression in HPGs and mandibular glands by increasing the bee ages.
Differential expression of CYP450 genes between two groups of the current study is due to the role of HPGs in the secretion of detoxification compounds during the nectar process. Although 9_DO HBs are considered nurse bees and are not in direct contact with environmental pollutants, but they have tasks such as cleaning the hive. Vannette et al. (Reference Vannette, Mohamed and Johnson2015) by comparing the expression of immune and detoxification genes in forager and nurse HBs showed that in parallel with increasing the age, the expression of many detoxification-related genes, particularly in the HPGs and mandibular gland are increased. Overall, it seems that the expression of genes encoding detoxification enzymes, signaling molecules involved in the JNK pathway and phagocytosis increased in parallel with aging.
Ribosome related genes
Active genes in ribosome biogenesis (WDR36, WDR43 and WDR3) have higher expression in 9_DO HBs and only 4 genes involved in this pathway (WDR 28 and WDR 3) had higher expressions in 1_DO HBs. It can be concluded from the results of current study that the higher expression of 39S ribosomal genes in 9_DO HBs is to meet the energy providing for increased cell activity by aging from day 1th to day 9th and controlling the cell growth rate that started from worker HBs formation. The higher expression of 40S and 60S ribosomal protein genes in 1_DO HBs revealed that the protein subunits make at the early stage of growth and by increasing the age, the produced subunits used for ribosomal biogenesis process. This can be the reason for lower and higher gene expressions of ribosomal protein and ribosome biogenesis, respectively, in 9_DO HBs. Due to the fact that the main cellular activities in adult workers include basal cell processes and the production of RJ, older HBs need to produce more protein for maintaining cells to increase number and size, and also producing RJ.
Proteasome related activity
The proteasome is responsible for the folding, sorting and degradation of proteins. Studies have shown that the ubiquitin proteasome system in social insects such as HBs is involved in aging-related mechanisms (Lee et al., Reference Lee, Lee and Min2015; Shih et al., Reference Shih, Huntsman, Flores and Snow2020). The higher activity of this process in 9_DO HBs may be due to the needs of 9_DO HBs to expand the range of their social tasks. It also is in line with higher protein production and cell maintenance activities in this group.
Vitamins related pathways
It is somewhat surprising that vitamins B5 and B3 have the highest vitamin share in RJ content and vitamin B6 has very low amount (Collazo et al., Reference Collazo, Carpena, Nuñez-Estevez, Otero, Simal-Gandara and Prieto2021). In the present study, the expression of genes in vitamin B6 metabolism had a higher expression in 9_DO HBs compared with 1_DO. This could be due to the feeding role of 9_DO HBs.
Aromatic amino acids
Phenylalanine and tyrosine are aromatic amino acids that their synthesis and degradation pathway are necessary for production downstream secondary metabolites such as defense compounds, pigment, flavonoids, neuroprotectants, vitamins, cofactors, and also hormones (Parthasarathy et al., Reference Parthasarathy, Cross, Dobson, Adams, Savka and Hudson2018). In recent research, the activity of Phenylalanine, tyrosine and tryptophan biosynthesis pathway (ame00400) was increased in 1_DO HBs compared with 9_DO. After eclosion (the emergence of a worker from the pupa), melanization and sclerotization occurs in cuticle of young workers. The most important event is the differentiation of the exocuticle that cause to harden the outer procuticle. Hardening (tanning/ sclerotization) is usually accompanied with darkening or melanization (Anderson, Reference Anderson, Kerkut and Gilbert1985; Hopkins and Kramer, Reference Hopkins and Kramer1992). Before sclerotization begin, the level of tyrosine increases in the hemolymph. Then it is converted to dihydroxyphenylalanine (Dopa) (Hopkins and Kramer, Reference Hopkins and Kramer1992; Hopkins et al., Reference Hopkins, Starkey, Xu, Merritt, Schaefer and Kramer1999) and it was also converted to N-acetyl-dopamine and transported via the pore channels to the epicuticle for quinones production (tanning agents). Ommochromes, derived from tryptophan, are important group of pigments that produce yellow, red, and brown colors after eclosion in insects (Gillott, Reference Gillott2005).
Neuroactive ligand-receptor interactions
Although the role of defective 3 homolog protein mentioned as ‘unknown function’ in annotation of HB genome, it enrolled in asymmetric cell division, polarized growth epithelial tight junction, cell-cell adherence junction and cell-polarity in various kinds of tissues (Manabe et al., Reference Manabe, Hirai, Imai, Nakanishi, Takai and Ohno2002; Wang et al., Reference Wang, Du, Fang, Yang, Zhang, Zeng, Ullrich, Lottspeich and Chen2006; Kunnev et al., Reference Kunnev, Ivanov and Ionov2009). Lower expression of Neuroactive ligand-receptor interaction genes is in agreement with those obtained by Han et al. (Reference Han, Fang, Feng, Hu, Hao, Ma, Huo, Meng, Zhang and Wu2017) that they examined the brain proteomes of newly emerged, worker, nurse, forager of Italian and the RJ bee strains. They reported the neuroactive ligand-receptor interaction protein of forager RJ bee on day 7th and 21th, nurse RJ bee on day 14th, Italian bee on day 7th and Italian forager bee on day 21th differentially up-regulated in comparison with other group. Given that the proteins resulting from the activity of this pathway not only produced in 7 to 21th days but also had a significant increase in its expressions, it is obvious that their transcription process has been done in previous stages. This concept is consistent with our finding that the expression of neuroactive ligand-receptor interaction pathway decreased from 1 to 9th day.
RNA degradation pathways
RNA degradation includes several pathways and a crucial role in mechanisms related to transcriptional regulation. It seems that by increasing intensity and number of vital activities in 9_DO HBs, they need a higher level of regulation of gene expression. HBs face various stresses including food stress resulting from changes in plant resources, temperature stresses due to their temperature adaptation to the environment inside the hive, and regulation of the temperature of the muscles during flight. Proteostasis is one of the most important processes that help HBs to adapt to such conditions via regulate protein synthesis, folding and degrading processes. The heat shock proteins have a key role in this complicated process. The Hsp10, Hsp90, Hsp83, Hsp70, Hsp40 and Hsp60 detected in HB genome. Among these, Hsp40 protein has a J-domain, and based on its conserve sequence, this protein divided into the DnaJA, DnaJB and DnaJC (Craig et al., Reference Craig, Huang, Aron and Andrew2006; Craig and Marszalek, Reference Craig and Marszalek2017).
The expression of heat shock protein genes in 9_DO HBs is lower than that of 1_DO and only Hsp60 had higher expression.
Based on the lower expression of the most of Hsp genes, it can be concluded that the both study groups did not experience any particular thermal stress. But the higher expression of Hsp40 in 9_DO HBs can be assumed as its function in protein folding, transport and degradation process in other pathways. Further studies regarding the role of Hsp40 in development of HBs in normal condition and without heat stress is strongly recommended.
PPI networks and hub genes
Results of comparing the head transcriptome of HBs at 1 and 9-daysold indicated that the increasing age of HB requires extensive changes in their gene expression profile as well as its PPI networks. The important issues as a general result of comparing transcriptome are the identification of key genes involved in certain characteristics, understanding the overall mechanism and identification of gene expression trends and for implying in future studies. In the current study, DEAD-box ATP-dependent RNA helicase 20, ATP-dependent RNA helicase DDX47, ATP-dependent RNA helicase DDX28, tRNA (cytidine(32)/guanosine(34)-2′-O)-methyltransferase and pseudouridylate synthase 7 have the most connections in the protein network, respectively. All five top identified genes have crucial roles in important cell process and also, are up-regulated in 9_DO HBs. DEAD-box ATP-dependent RNA helicase 20 gene codes the DEAD-box proteins. This protein is detected in all 3 life domains and carries out several roles in macromolecular machine assembly like the ribosome and spliceosome, controlling the gene expression, mRNA export from nucleus, self-splicing RNA introns folding (Jarmoskaite and Russell, Reference Jarmoskaite and Russell2011). Also, ATP-dependent RNA helicase DDX47 belongs to the DEAD-box family of ATP-dependent RNA helicases and has a role in apoptosis, rRNA transcription, mRNA splicing, pre-rRNA processing (Lee et al., Reference Lee, Rho and Chun2005; Sekiguchi et al., Reference Sekiguchi, Hayano, Yanagida, Takahashi and Nishimoto2006; Singleton et al., Reference Singleton, Dillingham and Wigley2007; Zhang et al., Reference Zhang, Forys, Miceli, Gwinn and Weber2011; Liu and Imai, Reference Liu and Imai2018). tRNA (cytidine(32)/guanosine(34)-2′-O)-methyltransferase acts as healthy growth, regulation of translation, RNA methylation via modification of tRNA the post-transcriptional modification of tRNA anticodon loop (Guy et al., Reference Guy, Podyma, Preston, Shaheen, Krivos, Limbach, Hopper and Phizicky2012). Pseudouridylate synthase 7 is catalyzer of RNA pseudouridylation (Safra et al., Reference Safra, Nir, Farouq, Slutskin and Schwartz2017; Guzzi et al., Reference Guzzi, Cieśla, Ngoc, Lang, Arora, Dimitriou, Pimková, Sommarin, Munita and Lubas2018; Shaheen et al., Reference Shaheen, Tasak, Maddirevula, Abdel-Salam, Sayed, Alazami, Al-Sheddi, Alobeid, Phizicky and Alkuraya2019). Interestingly, cytochrome P450 305a1 gene is one of 20 significant identified genes and plays a key role in the mechanisms of immunity, detoxification and removal of foreign compounds from nectar. Some of the selected genes are briefly discussed. In general, what stand out regarding in identified hub genes from PPI network is that these genes not only participate in four important activities of RNA structural modification, regulating gene expression and protein production, neutralizing potential stress and detoxification, but also their expressions increased by aging. In addition, it can be stated that the activity of HB head is focused on energy supply and increasing the protein production and export instead of storage, while actively trying to increase the level of cell dynamics and immunity for reducing stress side effects.
Conclusion
The transcriptome study of the HB's head during development can lead to further understanding about the growth, secretory capabilities of HPGs, RJ production, and identifying the key stage of certain behavior. This study set out to investigate the head transcriptome of 1 and 9 days old HBs. Results indicate that 9_DO nurse HBs up-regulate various basal pathways as amino acid and protein biosynthesis/metabolism, fatty acid biosynthesis/elongation/metabolism/degradation, starch and sucrose metabolism, monoscharid biosynthesis, oxidative phosphorylation, citrate cycle, glycolysis/gluconeogenesis, metabolic pathways, transporters to improve their development and preparing for doing their duties in hives. An important finding of our study is that the expression of MRJP genes has a higher expression in 1_DO HBs. By increasing the age and level of social activity, the necessity of 9_DO HBs increased for keeping the health of nurse HBs against the chemical compound of entered to the hive by foragers. Most of the process and attributes related to HPGs activities have an upward trend from 1_DO to 9_DO HBs and are age-dependent. The focus of 9_DO HB is on the fatty acid elongation and degradation process rather than fatty acid biosynthesis. The expression of aromatic amino acid and also vitamins genes are essential for early stage of life while. In 9_DO HBs, the energy supplying, reducing stress, protein production and export have a crucial role for support the body development and increase the social duties.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485322000554
Conflict of interest
The authors declare no competing interests
Data availability
The raw RNA-seq data have been deposited in the SRA database of database SRR14767741 and SRR14767740.
