Division of Labor among Worker Bees Is Associated with the Lipidomic Plasticity in Their Brains

Abstract

The division of labor is a dominant characteristic of honeybees and is accompanied by behavioral specialization and cognitive enhancement. As the central nervous system to control the labor-specific behaviors of honeybee, the brain is richest in lipid in terms of both diversity and abundance. In this study, an in-depth LC-MS/MS-based lipidomic method was applied to systematically characterize the brain lipid compositions of worker bees with three labor stages: newly emerged bee (NEB), nurse bee (NB), and forager bee (FB). A total number of 337 lipid species that assigned to 20 lipid classes were analyzed. The association of the brain lipidomes with the division of labors was suggested by the results of both the unsupervised and supervised multivariate pattern recognition analysis. More than 68% of the identified lipid species were found to be significantly changed in at least one comparison between NEB, NB, and FB. A total of 81 lipid species were identified as the potential labor-featured molecules with VIP > 1 and p-adj < 0.05. The labor-featured lipids of FA(18:2), FA(18:3), FA(26:0), PC(18:0_18:3), PS(18:1_18:1), SM(d38:1), CoQ10, and CoQ9, as well as their interactions with 12 behavior-related genes, including AmEST-6, AmFABP, AmE75, AmDGAT2, AmLSD1, AmNPC1, AmABCA1, AmNMDAR1, AmHTT, AmNOS, etc., were revealed by the further IPA analysis. These findings demonstrate for the first time that the brain lipidomes of worker bees are associated with the stable differences in their labors, which help understand the function of brain lipids on the labor-dependent behaviors of honeybees.

1. Introduction

The honeybee (Apis mellifera ligustica) is a typical social insect with an age-related division of labor or temporal polyethism [1,2]. Worker bees live on average 6 weeks during production seasons, during which their division of labor is based on their highly stereotyped behavioral maturation [3]. Although with less-developed brains, NEBs start by cleaning the cell and keeping the brood warm for about 3 days after their emergence. From the age 4 to 12 days, workers function as NBs to feed the larvae and care for the queen. Then, between the ages of 12 to 21 days, workers mainly engage in the processing and storing of food, the construction of the comb, the regulation of hive temperature, and the guarding of the entrance. During the late stages of their life cycle (day 21–45), worker bees switch from working in the hive to foraging outside for pollen, propolis, nectar, and water until they die. Honeybees with various laborious tasks displayed highly organized social behaviors that are associated with the activities of their brain, making it a perfect model to study the connection between brain function and social behaviors [1,4,5].

Worker bees’ experience inside and outside the colony and the evaluation of the temporal and spatial stimuli create memory traces in their brains [6], and the proboscis extension reflex (PER) assay was widely used to investigate the brain activities [7]. It has been found that the individual variations in the PER response thresholds of NBs, PFs (pollen foragers), and NFs (nectar foragers) to task-specific stimuli were linked with their behavior phenotypes [8]. It was also found that NEBs were very unresponsive to sucrose and displayed poor learning capabilities compared with adults who have relatively mature brain functions [9,10]. Another study using cGMP to make young bees forage found that their memory skills of remembering rewarding flowers were improved, suggesting that the brain functions of worker bees were determined by their behaviors [11].

The nervous system has the highest lipid content and complexity to support a number of key functions, including synaptogenesis, neurogenesis, insulation, etc. [12,13]. Glycerophospholipids, sphingolipids, and glycerolipids are three major lipid categories identified in the brain [14]. About half of the brain lipids were comprised of GPs, which are broadly categorized based on their head groups as well as the different combinations with fatty acyl chains. GPs not only provide stability, fluidity, and permeability of the cellular membrane but also interact with various membrane proteins and act as reservoirs for the second messengers for their proper functions [15]. SPs consisting of a sphingo base are ubiquitous structural components of membranes and are involved in cellular signaling, including the regulation of numerous ion pumps, channels, and transporters [16]. In both insects and mammals, memory formation and cognitive impairment were affected by the membrane fluidity that is associated with the degree of unsaturation of the membrane lipids [17,18]. Glycerolipids are carbohydrate-attached lipids that provide energy to all living things [19]. Moreover, cholesterol, although not present in insects, resides in the brain of mammals in the form of myelin [20].

Honeybee labor-specific behaviors such as nursing and foraging are composed of systematically organized social interaction and communication that are controlled by the central nervous system (CNS). Over the past two decades, it has been increasingly interesting to understand the regulatory mechanisms underlying the honeybee behaviors at the molecular levels [21,22,23]. A batch of genes was found to be involved in the labor-related behavior changes, for example, USP, NMDAR1, NOS, GR10, Egr-1, Hr38, and Kakusei [24,25,26,27,28]. A study aiming to identify the brain transcription factors (TFs) of forager behaviors found CWO, NF-kappa B, EGR, PAX6, and Hairy co-expressed with nearly half of the upregulated genes in FBs [28]. CWO and Hairy were recently verified to regulate forager behaviors based on the results of high-throughput behavioral tracking, brain gene expression, and chromatin accessibility profiles [29].

On the one hand, the lipid metabolism within the brain is tightly regulated to maintain the structure and function of the brain. On the other hand, lipids may signal environmental status to modulate brain activities [30]. In honeybees, the associations between behaviors or behavior-related genes with dietary lipids such as oleic acid, linlenic acid, and linolenic acid were identified. Oleic acid was found to trigger the hygienic behavior in honeybees by associating with OBP16 and OBP18 [31]; linolenic acid or linoleic acid depletion may cause foragers to dance faster to the pollen with “complementary” fatty acids to balance its nutritional requirement [32]. A recent study examined the lipidomic changes in the brain of 7-day-old bees after challenging with neurotoxic insecticide and found the lipid species of CL (18:3/18:1/14:0/22:6), TG (6:0/11:2/18:1), and eLPE (18:0e) linked with the self-grooming behavior [14].

Given the above, it has become increasingly apparent that lipid metabolism within the CNS contributes to the regulation of labor-related behaviors of honeybees. However, the collective characterization and quantification of pools of brain lipids are required to provide adequate information about the association between brain lipid plasticity and behavioral specialization. The aim of this study was to broadly characterize the brain lipid compositions of NEB, NB, and FB and provide new insights into the role of lipid metabolism in the regulation of social behaviors in honeybees. In this study, lipidomic profiles of the brains of NEB, NB, and FB were created using the UPLC-Q-Exactive Orbitrap MS/MS method, and the differences between these sample groups were investigated using both unsupervised and supervised multivariate pattern recognition methods. Then the candidate labor-featured lipid molecules, as well as their relative levels, were uploaded to Ingenuity pathway analysis (IPA) for network construction, by which the essential lipid molecules and their interactions with behavior-related regulatory genes were obtained. A total of 8 labor-featured lipids, including FA(18:2), FA(18:3), CoQ10, etc., as well as their interactions with 12 behavior-related genes, including AmEST-6, AmE75, AmNPC1, etc., were revealed by the further IPA analysis.

2. Materials and Methods

2.1. Sample Collection

Honeybee colonies used in this study were maintained at the apiary of the Institute of Apicultural Research, Chinese Academy of Agricultural Science, Beijing, China. The experiment was performed based on six colonies from the same queen bee as replicates with the same management. Frames with capped brood before eclosion were collected from the field colonies and incubated at 34 °C and 70% humidity. Newly emerged bees (NEBs) were collected when adult bees emerged from their cells within 24 h. The 10-day-old worker bees showed nurse behavior in comb cells collected as nurse bees (NBs). The foraging bees (FBs) were collected from the entrances of the hives, whose legs were loaded with pollen. The reason for choosing NEB, NB, and FB is that nursing and foraging are the most distinct behavioral specialization among social insects such as ants and bees [33], and NEBs were used as a control for they are commonly considered developmentally immature [34]. The collected bees were immediately frozen in liquid nitrogen and stored at −80 °C until the dissection of brain samples for the subsequent lipid extraction.

The individual bees were fixed on a beeswax dish with dry ice to keep cool. The antennas and head capsules were carefully removed using a pair of forceps to expose the entire brain. Then, the hypopharyngeal glands, the ocelli, as well as the compound eyes were removed. From each of the six colonies, five brains from the bees with the same age/labor were mixed as one sample. Thus for the three labor stages, a total number of 18 samples were used for lipids extraction and lipidomic analysis.

2.2. Lipids Extraction

Lipid extraction was performed immediately after the dissection of brain samples. Lipid samples were prepared according to the method of Folch et al. [35] with several modifications. Brain samples of the same labor group from each colony were mixed with 2 mL CHCl₃: MeOH (2:1, v/v) in glass tubes and fully homogenized for two minutes’ vortex on ice. The mixture was combined with 500 μL HPLC-water (1:4, v/v, water phase: organic phase) and homogenized by 25 s vortex and 2.5 min’ incubation on ice three times. After centrifugation at 3000 rpm for 15 min, the lower organic phase was collected and transferred into new glass tubes. The isolation was repeated one more time as above. Then the lipid samples were completely dried under a stream of nitrogen and stored in a −80 °C freezer as dry pellets.

2.3. UPLC-Q-Exactive MS/MS Methods for Lipids

The untargeted lipidomics profiling was performed on ultra-high-performance liquid chromatography (UPLC) Q-Exactive Orbitrap mass spectrometer (Thermo Fisher, Waltham, MA, USA) equipped with an electron spray ionization (ESI) probe. The lipid extracts were separated using a Cortecs C18 column (2.1 × 100 mm, Waters, Milford, MA, USA), and sample analysis was carried out in both negative and positive ionization modes. The column temperature was set at 40 °C, and the flow rate was 250 μL/min. The injection volume was 1 μL in positive mode and 2 μL in negative mode. The binary solvent system of mobile phase A was 10 mM ammonium acetate in HPLC-grade acetonitrile/HPLC-grade water (60:40 v/v), mobile phase B contained 10% HPLC-grade acetonitrile and 90% isopropanol (v/v). The samples were eluted with a linear gradient from 33% B to 98% B over 19.5 min, followed by 98% B for 10 min and re-equilibration with 33% B for 5 min.

The MS parameters for lipids analysis were as follows: the spray voltage was set at 3.2 kV in positive mode and 2.8 kV in negative mode. The capillary temperature is 320 °C, and the flow rate of sheath gas and aux gas is 35 arbs and 10 arbs, respectively. The data were acquired at mass ranges (m/z) 240–2000 for positive mode and 160–2000 for negative mode. The full scan and fragment spectra were collected at a resolution of 70,000 and 17,500, respectively [36].

To detect signal deviations related to instrumental drift and to assess the reproducibility of the data, quality control (QC) samples were prepared by pooling 10 μL of each of the NB samples and were inserted in an analytical batch after every six sample injections.

2.4. Lipid Identification

The lipids were identified and quantified using LipidSearch v4.1.30 (Thermo Fisher Scientific, San Jose, USA) based on the retention time and the characteristic fragment ions. A total of 5 ppm and 10 ppm mass tolerance was used for precursor and fragment, respectively. A retention time shift of 0.15 min was allowed for “alignment”. The lipids with an m-score less than 6 and chromatographic peak area >10⁷ were considered confident identifications. The peak area of lipids identified under the positive and negative modes was normalized separately by the total ion intensity. The normalized abundance of each lipid species was determined by dividing the peak area of each lipid by the sum of the peak intensity in that sample and multiplying by 10¹⁰. The formula was as follows:

Normalized abundance of lipid species x = [(Peak area of lipid species x)/(Peak area of all lipid species)] × 10¹⁰.

2.5. Statistical Analysis

The normalized abundance of all identified lipid species in positive and negative ion modes was imported into SMICA v14.1 (Umetrics, Umea, Sweden) for principal components analysis (PCA) and orthogonal projections to latent structures discriminate analysis (OPLS-DA). For both PCA and OPLS-DA, the scaling type was set as Pareto scaling, the scaling block was set as one, and the block weight was set as 1/sqrt. For OPLS-DA, a total of 200 permutations was used to estimate the over-fitting of the model, and the lipids with variable importance in projection (VIP) scores >1 were screened for further analysis.

For each of the comparisons, including NB_NEB, FB_NEB, and FB_NB, the significantly changed lipids were examined by unpaired two-tailed t-test (p-adjust < 0.05 by Holm–Sidak correction).

2.6. Unsaturation Parameters

The fatty acyl chains of GL classes, as well as the 5 dominant GP classes, including PE, PS, PC, PI, and CL, were used for the unsaturation analysis. The above GP classes are known to account for 99% of the total brain GPs of honeybees, thus denoted as “dominant” [14]. The ratio of the abundance of the saturated fatty acid (SFA), monounsaturated fatty acids (MUFA), as well as polyunsaturated fatty acids (PUFA) to the abundance of total fatty acids in GPs and GLs, were calculated as (with “*” representing “SFA”, “MUFA” or “PUFA”):

Ratio of * = (the abundance of */the abundance of total fatty acids) × 100%.

The unsaturation index (UI) of fatty acyl chains was calculated following Malaguti et al. [37], but using the relative abundance (RA) of the lipid species where the fatty acyl chains were from:

UI = 1 × (the RA of monoenoics) + 2 × (the RA of dienoics) + 3 × (the RA of trienoics)
+ 4 × (the RA of tetraenoics) + 5 × (the RA of pentaenoics) + 6 × (the RA of hexaenoics).

2.7. Ingenuity Pathway Analysis (IPA)

IPA tools (software version 52912811, QIAGEN Inc., Hilden, Germany) was used to explore the biological functions of the labor-featured lipid species. The lipids with the identification from the “human metabolome database” (https://hmdb.ca/, accessed on 1 March 2022) and their abundance were uploaded for network construction. The Fisher’s exact test was employed to calculate a p-value to determine the probability of the association between the lipids in the data set and the predefined relations with genes and endogenous chemicals in the Ingenuity Pathway Knowledge Base due to random chance alone.

2.8. qRT-PCR Analysis

TRIzol reagent (Invitrogen, Waltham, MA, USA) was used to extract the total RNA of the brain samples following the manufacturer’s protocol. The quantity and integrity of the RNA samples were examined by NanoDrop 2000 spectrophotometer (Thermo Fisher, Waltham, MA, USA) and 1.2% agarose gel. The cDNA was synthesized from a 1 μg RNA template using PrimeScrip RT reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, Kusatsu, Japan). A total of 0.8 μL diluted cDNA was then used for subsequent qPCR analysis in a final volume of 10 μL reaction system using TB Green Premix Ex Taq II (Tli RNaseH Plus) (Takara, Tokyo, Japan) (3 technological replicates). The reactions were conducted on LineGene 9600 Plus (Bioer, Hangzhou, China) using the following procedure: 95 °C for 30 s; 40 cycles of 95 °C for 5 s, and 58 °C for 30 s. The gene expression was calculated by the 2^−∆∆Ct method [38] using Ribosomal protein L32 (XM_006564315.3) as the reference gene. All the primer sequences were listed in Table S1.

3. Results

3.1. Lipidomic Plasticity in Brains of Worker Bees with Different Labors

Total ion current (TIC) chromatograms of the brain lipid extracted from NEB, NB, and FB are displayed in Figure S1. A total of 304 features in the positive mode and 240 features in the negative mode representing five categories were identified: Glycerophopholipids (GPs), Glycerolipids (GLs), Sphingolipids (SPs), fatty acyls, and Prenol lipids (PRs) (Figure 1a). These lipid features were assigned to 20 lipid classes, including Phosphatidylcholine (PC), Phosphatidylethanolamine (PE), Phosphatidylglycerolipid (PG), Phosphatidylinositol (PI), Phosphatidylserine (PS), Lyso-Phosphatidylcholine (LPC), Lyso-Phosphatidylethanolamine (LPE), Cardiolipin (CL), Phosphatidic acid (PA), Sphingomyelin (SM), Ceramide (Cer), Sphingosine (SPH), MonohexosylCeramide (Hex1Cer), DihexosylCeramide (Hex2Cer), Triacylglycerol (TG), Diacyglycerol (DG), Fatty Acid (FA), Acyl Carnitine (AcCa), Coenzyme Q (CoQ), and O-acyl-(ω-hydroxy)-fatty acid (OAHFA). Among them, PC, PE, CL, SM, TG, and FA were the most common and abundant lipid classes, with each of them composed of more than 20 species. GPs and SPs were detected in both ESI-positive mode and ESI-negative mode. In this study, the GP signals in the negative mode and SP signals in the positive mode were reserved for the quantitation. Because the fatty acyl chains of GPs were only able to separate in the negative mode and more stable instrument responses for SPs were observed in the positive mode [39,40]. Finally, a total of 337 lipid species were used for the subsequent quantification analysis (Tables S2 and S3).

To identify and characterize the possible associations between the division of labor and brain lipidomic plasticity, principal component analysis (PCA) was performed. As shown in Figure 1b, three clusters were clearly separated with PC1 (accounting for 49.5% of the total variance) and PC2 (accounting for 20.7% of the total variance). It was also shown that QC samples were clustered tightly in PCA scatter plots, indicating a stable performance of the UPLC-Q-Exactive MS/MS profiling.

The abundance variation of lipid classes among NEB, NB, and FB samples was investigated (Figure 1c–g). For the GP category, the abundance levels of PC, LPC, LPE, and CL were similar in NEB, NB, and FB samples. The level of PE was significantly higher in FB, while the level of PS showed an opposite trend. The abundance of PI continuously decreased from NEB to NB, then to FB. The abundance of PG was higher in NB than in the other two stages, while the abundance of PA was the lowest in NEB (Figure 1c and Table S4). For the SP category, SM, Hex1Cer, and Hex2Cer showed the lowest abundance in NEB. In contrast, Cer showed the lowest abundance in FB. The abundance of SPH was significantly higher in NB (Figure 1d). For the GL category, the abundance of both TG and DG was higher in NEB compared with NB and FB (Figure 1e). The abundance of free fatty acids was higher in both NB and FB (Figure 1f). The overall abundance of OAHFA and AcCa was not significantly different among the sample groups (Figure 1f). In addition, a limited number and low peak signals were detected for both of them (Tables S2 and S3). For the PR category, CoQ9 and CoQ10 were found with higher levels in FB compared with NEB and NB (Figure 1g).

3.2. The Degree of Unsaturation of Fatty Acyl Chains in GPs and GLs

The degree of unsaturation of fatty acyl chains is one of the most important features known to be vital for the fluidity of neural cell membranes [17,18]. In this study, the largest number of identified lipid species were assigned to the GP and GL category (Figure 1a), which constitutes the significant components of neural cells.

To investigate the degree of unsaturation of the fatty acid chains in GPs, the ratio of SFA, MUFA, PUFA to the total FA, and unsaturation index (UI) was calculated (Figure 2a,b). MUFAs accounted for the highest proportion of the total fatty acyl chains of GPs in all the brain samples, followed by SFAs and PUFAs (Figure 2a). The ratio of SFA to total FA were 30.44% ± 0.28%, 30.40% ± 0.32%, and 29.91% ± 0.37% in NEB, NB, and FB, respectively. The ratio of MUFA to total FA was 45.75% ± 0.82% in NEB, and significantly decreased to 40.34% ± 0.66% and 40.72% ± 0.65% in NB and FB, respectively (p < 0.05). The ratio of PUFA to total FA was 23.81% ± 0.68% in NEB, and significantly increased to 29.26% ± 0.50% and 29.36% ± 0.54% in NB and FB, respectively (p < 0.05). In addition, the UI of the fatty acyl chains in GPs was 1.08 ± 0.01 in NEB and significantly increased to 1.18 ± 0.01 in both NB and FB (p < 0.05, Figure 2b).

The percentages of the fatty acid residues with different numbers of double bonds (DBs) were further investigated: the fatty acid residues with more than two double bonds showed increased levels in NB (40.34% ± 1.21% in PEs, 41.58% ± 2.14% in PCs, and 40.51% ± 1.84 in PSs) and FB (40.56% ± 0.70% in PEs, 39.10% ± 1.53% in PCs, and 47.36% ± 2.07% in PSs) than in NEB (32.46% ± 1.34% in PEs, 28.50% ± 1.66% in PCs, and 34.88% ± 2.06% in PSs) (Figure 2c). It was noted that the fatty acyl chains with 18 carbons, including FA(18:0), FA(18:1), FA(18:2), and FA(18:3), accounted for the highest proportion of the total fatty acid residues: for PE, PS, and PC, the percentages were 87.05~87.89%, 99.62~99.67%, and 77.57~79.62%, respectively (Figure S2).

Differently from GPs, the ratio of SFA in the total fatty acid residues of GLs was 37.95% ± 2.16% in NEB and greatly decreased in NB (18.79% ± 3.30%) and FB (20.20% ± 1.99%); the proportion of MUFA and PUFA displayed a reversed wave pattern, with the lowest level of MUFA (34.45% ± 3.51%) and the highest level of PUFA (46.76% ± 6.75%) appeared in NB (Figure 2d). The index UI of the fatty acid residues in GLs displayed a similar pattern to PUFA (Figure 2e). The greatest differences were seen for the proportion of PUFAs, which increased about 28% from NEB (18.29% ± 2.27%) to NB (46.76% ± 6.75%) (p < 0.05) (Figure 2d).

The distribution of the number of DBs of the fatty acid residues of DGs and TGs is shown in Figure 2f. In DGs, the fatty acid residues with two DBs account for the highest percentage in all labor groups (57.70% ± 2.28% in NEB, 34.27% ± 4.11% in NB, and 54.37% ± 5.00% in FB). However, a more diverse distribution of DBs in the fatty acid residues of TGs was observed, with the number ranging from 0 to 9. In both TGs and DGs, the dominant fatty acid residues are 18Cs with DB numbers ranging from 0 to 3 and 16Cs with DB numbers of 0 and 1. In TGs, the saturated fatty acid residues with 12C and 14C were also detected with relatively strong signals (Figure S3).

3.3. The Compositions of Free FA in the Brain of Worker Bees with Different Labors

Free fatty acids, particularly polyunsaturated fatty acids (PUFAs), play critical roles in the regulation of both the structure and the function of neurons, glial cells, and endothelial cells [41]. There are 16 saturated fatty acids (SFAs), 6 monounsaturated fatty acids (MUFAs), and 6 PUFAs identified in the brain of worker bees. SFAs with an even number of carbons from 16 to 30 were identified to be dominant saturated fatty acids in the brain of worker bees (Figure 3a). A total of 6 SFAs with an odd number of carbons identified with pretty weak signals, including FA(29:0), FA(27:0), FA(25:0), FA(23:0), FA(17:0) and FA(15:0) (Figure 3a). These observations are consistent with the findings that FAs with an odd number of carbons are uncommon in animals [42].

The ratio of MUFA was higher (23.1% ± 2.0%) in NEB and close in NB (18.15% ± 3.55%) and FB (16.95% ± 2.51%), which is consistent with the findings in GP (Figure S4 and Figure 2a). The most abundant MUFA was FA(18:1), which is a lot more than the others, followed by FA(16:1) and FA(20:1) (Figure 3b). For all the free FAs in the brain, the percentage of PUFAs was 7.43% ± 0.54% in NEB, 19.35% ± 3.04% in NB, and 13.82% ± 2.35% in FB (Figure S4). The dominant PUFAs were FA(18:2) and FA(18:3), which accounted for more than 92% of the total PUFAs. The level of both FA(18:2) and FA(18:3) increased more than two folds in NB and FB compared with NEB (Figure 3c). In addition, the accumulation of PUFAs in NB and FB was identified for 5 out of the 6 PUFAs, with the exception of FA(20:3) (Figure 3c).

3.4. Identification of Labor-Featured Lipid Species

The lipidomic data of the brain samples from NEB, NB, and FB were analyzed using both unsupervised and supervised multivariate pattern recognition methods. At the exploratory stage, principal component analysis (PCA) was used and suggested a very clear group clustering based on the division of labor (Figure 1b). Then multivariate statistical analysis was performed using orthogonal partial least-squares discriminative analysis (OPLS-DA) (Figure S5a–c). In our predictive models, all the values of R²X (cum), R²Y (cum), and Q² (cum) were larger than 0.5, suggesting suitable prediction abilities and classifications to distinguish the samples. The subsequent permutation tests (n = 200) further validated the robustness of all the models (Figure S5d–f). Candidate lipids associated with different labors within each of the comparisons were screened based on the VIP score from the OPLS-DA model (Figure S5g–i) together with the adjusted p-value from Student’s t-test (VIP > 1, p-adj < 0.05).

We found 51 candidate lipids between NB and NEB, 66 candidate lipids between FB and NEB, and 16 candidate lipids between FB and NB (Figure 4 and Table S5). As shown in Figure 4, a total of 40 lipids were found to simultaneously either increase or decrease in both NB and FB compared with NEB, an indication of the probability of their association with labor plasticity or complexity. Interestingly, PE(18:3_18:3) and SM(d38:2) were found universally changed in all of the three comparisons, with the level of PE(18:3_18:3) increasing 2.9-fold from NEB to NB and 3.7-fold from NEB to FB. The level of SM(d38:2) displayed a similar dramatically upward trend with the fold-change number of 3.5 and 5.6 in the comparison of NB_NEB and FB_NEB, respectively. The lipids that were always positively (+) or negatively (−) associated with a certain labor were listed. It was also noted that 8, 17, and 6 lipid species identified as featured lipids uniquely for the comparisons of NB_NEB, FB_NEB, and FB_NB, respectively. The average abundance, fold-change number, and significance text for each of the lipids were provided in Table S6.

3.5. Genes Associated with the Labor-Featured Lipid Species

Ingenuity pathway analysis (IPA, QIAGEN Inc.) is a powerful web-based tool for systematic analysis of genes and metabolites involved in the biological events by identifying the enriched canonical networks, thus is productive for biomarker screening. To explore the possible roles of the candidate lipids in the labor-featured behaviors, their HMDB identifiers (mapping rate > 76%) and the relative abundance were uploaded to IPA (input data referred to in Table S7) to extract the interactions with genes in mammals. Their homologous genes in honeybees (Apis mellifera, Taxid: 7460) were obtained by protein sequence blast as well as gene function annotation in the NCBI database (Table S8).

The first networks of labor-featured lipid species in NB_NEB, including PC(18:0_18:3) (18:0/18:3 phosphatidylcholine), FA(26:0) (cerotic acid), FA(18:2) (linoleic acid) and FA(18:3) (linolenic acid), showed connected with 4 behavior-related genes, including AmEST-6 (homologous to CES1), AmDGAT2 (homologous to DGAT2), AmFABP (homologous to FABP) and AmE75 (homologous to PPARG) (Figure 5a). FA(18:2) has an “activation” relationship with AmEST-6; “expression” relationship with AmDGAT2; “chemical-protein interactions”, “expression”, and “transcription” relationship with AmFABP; and “chemical-protein interactions”, “expression”, and “activation” relationships with AmE75. FA(18:3) has “chemical-protein interactions”, “expression”, and “translocation” relationship with AmFABP; and “chemical-protein interactions” with AmE75 (Table S9).

Two networks were obtained from the comparison of FB_NEB (Figure 5b,c). The first network containing PC(18:0_18:3), FA(18:2), FA(18:3), and 1,2-dioleoylphosphatidylserine (PS(18:1_18:1)). The importance of AmE75 was emphasized again for its connections with more than 60 genes or chemicals in this network. There are seven behavior-related genes in this network, including AmEST-6, AmFABP, AmDGAT2, AmE75, and AmLSD1 (homologous to PLIN2) (Figure 5b and Table S9). AmLSD1 has an “expression” relationship with FA(18:2). The relations between the featured lipids PC(18:0_18:3), FA(18:2), and FA(18:3) with the other four genes as described above. The second network from IPA results was associated with C20 sphingomyelin (SM(d38:1)). Both AmNPC1 (homologous to NPC1) and AmABCA7 (homologous to ABCA1) have “expression” relations with SM(d38:1) (Figure 5c and Table S9).

The lipid Coenzyme Q10 (CoQ10), which has the greatest VIP value in the FB_NB group, was highlighted in the fourth network (Figure 5d). CoQ10 is connected with 18 genes or chemicals that are involved in behaviors of cognition, learning, and memory. Among the above genes, eight genes have been found to be involved in learning and memory, including AmNMDAR1 (homologous to GRIN1), AmNMDAR2 (homologous to GRIN2A), AmHTT (homologous to HTT), AmBUFFY (homologous to BAX/BCL2), AmNOS (homologous to NOS2/NOS3), and AmE75. These genes all have “expression” relationships with CoQ10. In addition, AmHTT also has an “expression” relationship with CoQ9 (Ubiquinone 9).

The gene expression changes of the candidate genes associated with labor-featured lipids were validated using qRT-PCR of the same brain samples in NEB, NB, and FB. Of the 12 candidate genes, 10 genes showed significantly differentially expressed in at least one comparison of NB_NEB, FB_NEB, or FB_NB (Figure 6). The expression of AmEST-6 and AmE75 significantly (NB_NEB) increased by 2.0-fold (p = 0.02) and 2.4-fold (p = 0.000003) from NEB to NB, respectively (Figure 6a). While the expression of AmFABP significantly (p = 0.00006) decreased by 39% from NEB to NB. AmDGAT2, AmLSD1, AmE75, AmNPC1 and AmABCA1 were significantly (p < 0.03) upregulated from NEB to FB, and the fold-change varied from 1.2 to 5.9 (Figure 6b). The expression of AmNMDAR1, AmHTT, and AmNOS were significantly (p < 0.02) increased by 1.7-, 1.3-, and 1.4-fold from NB to FB (Figure 6c). It is notable that AmE75 showed significantly upregulated by 2.4-fold (p = 0.000003), 5.9-fold (p = 0.000005), and 1.9-fold (p = 0.0002) in the comparisons of NB_NEB, FB_NEB, and FB_NB, respectively.

4. Discussion

In honeybees, the brain is closely associated with learning, memory, and complex social behaviors for its function in processing second- or higher-order inputs from all sensory organs [21]. The brain is also known as the richest organ in terms of the content and diversity of lipids [12]. Many studies in mammals revealed that brain lipid compositions correlated with memory functions and cognitive development [43,44,45]. In honeybees, it has been found that certain lipid classes or species were enriched in the brain more than in other organs [14,46].

In the present study, adult worker bees with three different labor stages were used: (1) NEBs (developmental immature); (2) NBs (take nursing jobs in-hive); (3) FBs (foraging outside the hive). Lipidomic plasticity of NEBs, NBs, and FBs were monitored via UPLC-Q-Exactive MS/MS technologies, and a total number of 20 lipid classes assigned to 5 major lipid categories were comprehensively characterized (Figure 1a and Tables S2 and S3).

4.1. Dynamic Changes of Lipids in the Brain

The significant changes in the dominant membrane GPs and SPs across the transitions among NEB, NB, and FB were seen multiple times at the class levels (Figure 1c,d). This diversity may be interpreted by the investigation of more than 60% of the lipid metabolism-related genes were found differentially expressed in the brain of worker bees with different labors, including the key enzymes catalyzing the de novo biosynthesis of PE (CK/EK, PCYT, CPT/EPT), the switch between PC to PS (PSS1) and PS to PE (PSD) [23]. Our lipidomic results revealed that the level of PE was significantly increased in FB than in NB (p < 0.05), and the up-regulation of PSS1 and PSD in the brain of FB than in NB were verified by qRT-PCR (data not shown). It is also noted that the variance of PE, PS, PC, and PI among different labors was unlikely to the whole body analysis by Martin et al. [47], an indication of the spatio-temporal specificity of lipid composition in honeybees.

4.2. Unsaturation of Fatty Acids in the Brain

The unsaturation of FAs, either embedded in membrane lipids or exist in a free state, is an important index related to lipid peroxidation in various organisms, including honeybees [48]. Here we found that the fatty acyl chains from GP, GL, and free FA had the increased levels of unsaturation in NB and FB compared with NEB (Figure 2, Figure 3, and Figures S2–S4), providing new ideas for understanding the link between behavior diversity and the degree of the sensitivity to oxidation that determined by the DBs in the fatty acyl chains [47,48]. The distribution of chains with different numbers of carbons and DBs revealed that the increased unsaturation was greatly contributed by linoleic acid (FA(18:2)), an omega-6 FA, and linolenic acid (FA(18:3)), an omega-3 FA (Figures S2, S3 and Figure 3). Honeybee brain has a much higher proportion of FA(18:3) to the total FA than the other body parts, which is the opposite of that of FA (18:2) [46]. Arien et al. also found that the positive correlation between diet FA(18:3) and brain FA (18:3) was slight but significant (p = 0.028), and FA(18:3)-poor diets resulted in the reduced learning abilities and reduced size of the hypopharyngeal gland (an important nurse organ) of adult bees [46]. According to our lipidomic results, the highest levels of both FA(18:3) and FA(18:2) were detected in NB, followed by the second higher levels detected in FB, indicating its association with the nurse behavior and complicated brain activities (Figure 3). The above results were consistent with the findings of FA(18:3) and FA(18:2) as the potential labor-featured lipid species for their associations with behavior-related genes (Figure 5).

4.3. Brain Lipids Associate with Labor-Featured Behaviors

In the present study, a total of 81 candidates featured lipids were obtained based on the OPLS-DA (VIP > 1) and t-test (p-adj < 0.05) (Table S5). In an earlier study, Morfin et al. found a positive correlation between CL(18:3/18:1/14:0/22:6), TG(6:0/11:2/18:1), and LPE 18:0e and the behavior of intense self-grooming [14]. Thus, we focus on the lipids that are stably and positively linked with specialized tasks. Including DG(18:3_18:3), for it was found upregulated in every comparison that NB involved (VIP = 2.96 for FB_NB, and 2.16 for NB_NEB), thus considered a biomarker for NB behaviors. Co(Q10), PC(18:0e_18:1), PE(18:1_18:2) and SM(d36:2), with VIP ranges from 1.08 to 5.03, were suggested as the possible markers for FB behaviors (Table S6). In addition, 12 TGs, 4 DGs, 4 Pes, and 2 PCs, which were positively associated with NEB, were more important for early brain development than the complex information processing required by both in-hive and out-of-hive jobs (Table S6). More broadly, PE(18:3:18:3) and SM(d38:2) significantly varied between any two of the labor groups with a value of VIP > 1.31, implying the probability of their association with labor plasticity or complexity (Table S6).

4.4. Brain Lipids Associate with Behavior-Related Genes

According to IPA networks constructed based on the measurement of labor-featured lipids, 25 behavior-related genes, with 60% of them having homologs in honeybees, were found to directly or indirectly interact with the candidate labor-featured lipids (

Table 1. The relationship between behavior-related genes and labor-featured lipids. “-“ indicates “data not available”.

Lipid Species	Interactive Network	Gene Symbol	Behavior	Relationship (Lipid to Gene)	Relationship (Gene to Lipid)	Gene Symbol in Honeybees a
CoQ(10)	FB_NB	PPARA	feeding	expression	-	AmE75
CoQ(10)	FB_NB	NOS2	-	expression	-	AmNOS
CoQ(10)	FB_NB	NOS3	cognition, learning, and memory	expression	-	AmNOS
CoQ(10)	FB_NB	GRIN1	cognition, learning, and memory	expression	-	AmNMDAR1
CoQ(10)	FB_NB	GRIN2A	cognition, learning, and memory	expression	-	AmNMDAR2
CoQ(10)	FB_NB	HTT	cognition, learning, and memory	-	expression	AmHTT
CoQ(10)	FB_NB	BAX	cognition, learning, and memory	expression	-	AmBUFFY
CoQ(10)	FB_NB	BCL	learning and memory	expression	-	AmBUFFY
CoQ(10)	FB_NB	IL1B	cognition, learning, and memory	expression	-	-
CoQ(10)	FB_NB	Pro-inflammatory Cytokine	cognition	localization	-	-
CoQ(10)	FB_NB	SNCA	cognition, learning, and memory	localization	-	-
CoQ(10)	FB_NB	IL6	cognition, learning, and memory	expression	-	-
CoQ(10)	FB_NB	AGER	learning	expression	-	-
CoQ(10)	FB_NB	INS	cognition, learning, and memory	localization	-	-
CoQ(9)	FB_NB	HTT	cognition, learning, and memory	-	expression	AmHTT
CoQ(9)	FB_NB	IL1B	cognition, learning, and memory	expression	-	-
CoQ(9)	FB_NB	IL6	cognition, learning, and memory	expression; localization	-	-
FA(18:2)	NB_NEB; FB_NEB	CES1	feeding	activation	-	AmEST-6
FA(18:2)	NB_NEB; FB_NEB	DGAT2	feeding	-	expression	AmDGAT2
FA(18:2)	NB_NEB; FB_NEB	PPARG	feeding	chemical-protein interactions; activation	chemical-protein interactions; expression	AmE75
FA(18:2)	NB_NEB; FB_NEB	FABP	feeding	chemical-protein interactions; expression; transcription	-	AmFABP
FA(18:2)	NB_NEB; FB_NEB	FADS2	feeding	molecular cleavage	-	-
FA(18:2)	NB_NEB; FB_NEB	FFAR1	feeding, cognitioning, linoleic acid preference, oleic acid preference	activation	-	-
FA(18:2)	FB_NEB	PLIN2	feeding	expression	-	AmLSD1
FA(18:3)	NB_NEB; FB_NEB	FADS2	feeding	molecular cleavage	-	-
FA(18:3)	NB_NEB; FB_NEB	PPARG	feeding	-	chemical-protein interactions	AmE75
FA(18:3)	NB_NEB; FB_NEB	FABP	feeding	chemical-protein interactions; expression; transcription	-	AmFABP
FA(18:3)	NB_NEB	FFAR4	feeding, cognition, linoleic acid preference	activation; phosphorylation	-	-
PC(18:0_18:3)	NB_NEB; FB_NEB	PPARG	feeding	-	expression	AmE75
SM(38:1)	FB_NEB	NPC1	walking	-	expression	AmNPC1
SM(38:1)	FB_NEB	ABCA7	spatial learning	-	expression	AmABCA1

^a The homologous genes in honeybees were obtained from protein blast (https://blast.ncbi.nlm.nih.gov/Blast.cgi, accessed on 3 May 2022) and the annotation in the published papers.

). A percentage of 83% of the honeybee homologs were identified to be differentially expressed between the labor groups (Figure 6). Based on the relationships between lipids and genes (Table 1), the above genes could be classified into three groups:

(1)

The genes activated by labor-featured lipids. Both AmEST-6 and AmE75 were predicted to be activated by FA(18:2) (Table 1), and both were involved in juvenile hormone (JH) signaling. JH has long been known to regulate behavioral maturation and plasticity in insects [49,50]. EST-6 is a juvenile hormone esterase. In ants, it was known to be highly expressed in the brain of a nurse [51], and the oscillation in est-6 expression correlated with the labor shifting from nurse to forager [33]. In our study, the expression of AmEST-6 was upregulated in NB than in NEB (p = 0.02), while no significant difference was seen between FB and NEB (p = 0.07). Suggesting that AmEST-6 was more related to the nursing behaviors in honeybees.

Insect E75 is a nuclear hormone receptor especially important for JH responses during development [52]. E75/USP (ultraspiracle) is the counterpart of mammal PPAR/RXR (retinoid X receptor) complex [53,54,55], which plays roles in fatty acid transport and metabolism [56]. In honeybees, the transcription factor USP influences behavioral maturation, including the transition of in-hive labor to foraging outside and the expression of a batch of genes related to JH signaling, including E75 [55]. It was also found that RNAi of E75 resulted in the loss of locomotor rhythms in both firebrat and drosophila [57,58. In our study, the expression of AmE75 was correlated with the levels of labor-featured lipids FA(18:2), FA(18:3), PC(18:0_18:3), as well as Co(Q10), and suggested interacting with them through “activation”, “expression”, and “chemical-protein interaction” (Table 1). The above results indicated that E75 functions as effector genes of lipids for the regulation of insect behaviors.

(2)

The genes affected by labor-featured lipids at expression. The expression of 4 behavior-related genes, including AmNMDAR1, AmNOS, AmFABP, and AmLSD1, were regulated by labor-featured lipids and differentially expressed in the brain of worker bees with different labors. The expression of both NMDAR1 and NOS was affected by Co(Q10) in mammals. NMDR1 encodes the critical subunit of N-methyl-D-aspartate receptors, which plays a key role in the plasticity of synapses [25]. NOS is nitric oxide synthase that catalyzes the production of the signaling molecule nitric oxide (NO) from L-arginine [26]. Both NMDAR1 and NOS have been found to relate to learning and memory behaviors in honeybees [25,26]. However, the mechanism is not known yet. We here propose a possible role of Co(Q10) in brain function by affecting NMDAR1 and NOS. The gene FABP and AmLSD1 were possibly regulated by FA(18:2) and FA(18:3) (Table 1). FABP is a critical fatty acid-binding protein, was known to bind to FA(18:2) and FA(18:3) and is involved in the regulation of brain activity such as sleeping time and long-term memory in drosophila [59,60]. The association between FABP and FA(18:2), FA(18:3) was observed in this study, indicating the role of lipid FA(18:2) and FA(18:3) in the brain functions of honeybees.

(3)

The genes that affect the expression of labor-featured lipids. These gene-lipid interactions were predicted for E75, NPC1, HTT, ABCA1, DGAT2 and their downstream lipid effectors including FA (18:2), FA(18:3), PC(18:0:18:3), SM(38:1), Co(Q9) and Co(Q10). ABCA1 and DGAT2 have not been found to relate to behaviors in honeybees. However, they were known to be involved in spatial learning and feeding in mammals, implying their counterpart in honeybees may have similar roles in brain functions. Previous studies revealed that the metabolism of SM is associated with behavioral maturation by affecting the phosphorylation of CREB1 [28,29,61]. Here we found that an SM lipid—SM(d38:1)—might be regulated by a transporter for cholesterol and other lipids named NPC1, whose deficiency results in motor impairment, decreased exploratory activity, and cognitive impairment in both mice and drosophila [62,63]. Considering that the expression of NPC1 was upregulated in FB than in NB, we suppose that SM(d38:1) is a featured lipid for forager-related movement and cognition behaviors as the result of interacting with NPC1. Huntingtin (HTT), known to associate with the defensive behavior of worker bees [64], is a ubiquitously expressed nuclear protein that binds to a number of transcription factors, thus playing various roles in nerve cells [65]. The identification of HTT effected on CoQs in mammals gives us an idea of the lipids that may be involved in numerous brain activities as the downstream signals of HTT.

5. Conclusions

This is the first in-depth lipidomic profiling of the brain of worker bees that newly emerged with less-developed brain functions (NEBs), bees with nursing behaviors (NBs), and foraging behaviors (FBs). For the identification of possible behavior-featured lipids as well as their associations with brain gene expression, multiple methods, including PCA, OPLS-DA, IPA, and qPCR, were applied. A total of eight lipid species, including FA(18:2), FA(18:3), FA(26:0), PC(18:0_18:3), PS(18:1_18:1), SM(d38:1), CoQ9, and CoQ10 were implied to be promising indicators of labor-featured behaviors. The present study provides a lipidomic insight into the division of labor and behavioral specialization of honeybees; further experiments are required to determine their role in behavior regulation based on their interaction with proteins coded by Est-6, E75, FABP, NPC1, HTT, PI3K, etc.