Abstract

   Apis mellifera and Apis cerana are two important honey bee species
   widely kept and studied. But their queens differ greatly in egg-laying
   capacity. To determine the mechanisms of this difference, we compare
   gene expression, chromatin accessibility and spatial localization of
   differential genes in the ovaries of the two species in virgin queens
   and laying queens using ATAC-seq, RNA-Seq, homologous gene alignment
   and spatial transcriptome. The results reconfirm that the egg-laying
   capability of A. mellifera queens is significantly higher than that of
   A. cerana queens. The chromatin accessibility and nutrient cells ratio
   of A. mellifera queens are higher than those of A. cerana queens.
   Further investigations reveal that agxt2l (LOC408817) is significantly
   over-expressed in the ovaries of A. mellifera queens compared to A.
   cerana queens and is crucial for ovary development. Moreover, agxt2l
   can increase the phospholipid content in ovarian nutrient cells through
   the glycerophospholipid metabolism pathway to promote embryo formation
   and is regulated by brc-z1 (LOC552255). These findings suggest that the
   brc-z1-agxt2l signal pathway causes increased egg-laying in the queens
   of A. mellifera compared to the queens of A. cerana by accelerating
   lipid synthesis due to heightened glycerophospholipid metabolism.

   Subject terms: RNAi, Cell lineage, Reporter genes
     __________________________________________________________________

   The brc-z1-agxt2l signal pathway causes increased egg-laying in the
   queens of A. mellifera compared to the queens of A. cerana by
   accelerating lipid synthesis/transport due to heightened
   glycerophospholipid metabolism.

Introduction

   Honey bees (genus Apis) are highly eusocial insects that have benefited
   humans for centuries^[36]1. Beekeeping, the culture of both Apis
   mellifera and Apis cerana, has been used widely to produce bee
   products^[37]2–[38]4 and provide crop pollination^[39]5. Honey bees
   also help protect biodiversity^[40]6 and benefit human
   health^[41]7,[42]8.

   The biology of A. mellifera has been well studied^[43]9,[44]10.
   Compared to A. mellifera, the biology of A. cerana in Asia is less
   studied^[45]11 and many questions remain. There are many differences
   between the two sister species of honey bees. A. cerana has smaller
   body sizes than A. mellifera, but is better at collecting sporadic
   nectar in mountainous or forested areas^[46]12. A. cerana also has
   natural resistance against Varroa destructor^[47]13, better ability to
   cope with extreme climate^[48]14, and the ability to fight against
   giant hornets such as Vespa mandarinia^[49]15. A. cerana does not
   collect propolis^[50]16 and also fan in their nest entrances in the
   opposite direction, toward the hive^[51]17 instead of drawing air out
   by A. mellifera. A. mellifera has the advantage of higher honey
   production, which is the main reason for its spread throughout the
   world including its introduction to China^[52]18. However, the largest
   difference, and one that impacts honey yield drastically, is that A.
   cerana has much smaller colonies^[53]19. This is partly because they
   swarm more often^[54]20, partly because their queen egg-laying
   capability is much lower, compared to that of A. mellifera.

   The difference in queen egg-laying capacity between the two species is
   about 4-5 times. A. mellifera queens can lay 1500 to 2000 eggs per
   day^[55]21–[56]24, while A. cerana queens lay only 420 to 600 eggs per
   day^[57]25–[58]27. There are several possible mechanisms to cause this
   huge difference. One is that A. cerana queens have fewer ovarioles in
   their ovaries, with only 80-85% of the number in A. mellifera queens.
   However, the difference in ovariole count (across both ovaries) alone
   cannot fully explain the disparity in egg-laying capacity, because the
   ratio of ovarioles in queens of A. mellifera to A. cerana is only about
   1.5 to 2.0: A. mellifera queens show an average of 300 ovarioles^[59]28
   or a maximum of 400^[60]29, while A. cerana queens average of
   200^[61]30. The possibility that A. mellifera queens are fed much more
   often and/or with more or better food compared to A. cerana queens is
   not investigated. While such data do exist for A. mellifera
   queens^[62]31,[63]32, comparative data are missing in A. cerana. But
   even if the feeding rates are different, it still would require each
   ovariole to be more efficient, such that egg production per ovariole
   would be at least 2.5 times higher in A. mellifera than that in A.
   cerana.

   We hypothesize that A. mellifera queen ovarioles are able to produce
   mature eggs faster than those of A. cerana. This would require
   differences in physiology, which can be tracked via differences in gene
   expression. In this study, we used ATAC-seq, RNA-seq, spatial
   transcriptomics, and cross-species sequence comparisons to compare the
   gene expression differences in the ovaries of virgin queens and laying
   queens between the two honey bee species. We identified agxt2l as a
   potential gene: its higher expression in trophocyte cells enhances
   glycerophospholipid metabolism, thereby accounting for the higher queen
   egg-laying rate in Apis mellifera compared to A. cerana.

Results

Difference in queen egg-laying power between A. mellifera and A. cerana

   The number of eggs laid daily by A. mellifera queens was significantly
   higher than those in A. cerana (N = 10, Z = − 3.82, P < 0.01)
   (Fig. [64]1c). Body weight (N = 10, t = 4.32, df = 11.39, P < 0.01)
   (Fig. [65]1d) and the number of ovarioles (N = 10, t = 16.14, df = 18,
   P < 0.01) (Fig. [66]1e) in A. mellifera queens were significantly
   higher than those in A. cerana. The number of eggs laid per milligram
   of body weight per day (N = 10, t = 12.11, df = 18, P < 0.01)
   (Fig. [67]1f) and the number of eggs laid per ovariole per day in A.
   mellifera were significantly higher than those in A. cerana (N = 10,
   t = 20.80, df = 18, P < 0.01) (Fig. [68]1g). Details of statistical
   analyses are in Supplementary Data [69]1.

Fig. 1. Different egg-laying performance of queens in A. mellifera and A.
cerana.

   [70]Fig. 1
   [71]Open in a new tab

   a The newly laying queen in A. mellifera colony. b The newly laying
   queen in A. cerana colony. c The eggs laid daily of the newly laying
   queens of A. mellifera and A. cerana. d The body weight of the newly
   laying queens of A. mellifera and A. cerana. e The number of ovarioles
   in one ovary of A. mellifera and A. cerana. f The number of eggs laid
   per day per milligram of body weight. g The amount of eggs laid per
   ovariole per day in A. mellifera and A. cerana.

Sequence data quality control

   In this study, ATAC-seq and RNA-Seq analyses were conducted for the
   virgin queens (A. mellifera, AmV; A. cerana, AcV) and the newly laying
   queens (about two months old) (A. mellifera, AmQ; A. cerana, AcQ).
   After quality control, each sample obtained an average of 6.52 Gb of
   clean data (43.48 million clean sequence reads) (Supplementary
   Data [72]2). The Q30 values of each sample in the two omics were
   greater than 90% (Supplementary Data [73]2). For spatial transcriptome,
   2046 spot amounts were detected in the ovaries of A. mellifera queen
   (AmQ). For A. cerana queens (AcQ), 1625 spots were detected in the
   ovaries (Supplementary Data [74]2). These results showed that sequence
   quality was acceptable and the bioreplication was reliable.

DEGs between ovaries of A. mellifera and A. cerana queens

   DEGs were identified based on two sets of reference genomes (the
   reference genomes of A. mellifera and A. cerana). Pearson correlation
   coefficients were > 0.95 for all A. mellifera bioreplicates using the
   A. mellifera reference genome (Supplementary Fig. [75]1a), and > 0.91
   for all A. cerana bioreplicates using the A. cerana reference genome
   (Supplementary Fig. [76]1b).

   Because the reference genomes were different between A. mellifera and
   A. cerana, the expression level of the same gene may be different, even
   if they were the same, if compared using different reference genomes.
   To avoid this problem, genes showing the same expression trend in a
   comparison using two reference genomes were selected. In the comparison
   AmV vs AcV, there were 8481 and 10,315 DEGs using the A. mellifera and
   the A. cerana reference genome, respectively (Fig. [77]2a, d). Among
   these genes, 5976 genes were homologs (Fig. [78]2b). Among these
   homologs, 1154 DEGs showed the same expression trend under the dual
   species genomes and most of them were up-regulated in the virgin queens
   (A. mellifera, Fig. [79]2e). In the comparison AmQ vs AcQ, there were
   7341 and 9468 DEGs based on these two reference genomes respectively
   (Fig. [80]2a, d). Among these genes, 4821 genes were homologs
   (Fig. [81]2c). Among these homologs, 1543 DEGs showed the same
   expression trend in both bee species and most of them were up-regulated
   in the newly laying queens (A. cerana, Fig. [82]2f).

Fig. 2. Transcriptome data analysis results.

   [83]Fig. 2
   [84]Open in a new tab

   a The numbers of DEGs and the up-down-regulation relationship in the
   two comparisons (AmV vs AcV, AmQ vs AcQ) using the A. mellifera
   reference genome. b The Venn diagram of overlapped DEGs between AmV vs
   AcV (using the A. mellifera reference genome, AmV-m vs AcV-m) and AmV
   vs AcV (using the A. cerana reference genome, AmV-c vs AcV-c). c The
   Venn diagram of overlapped DEGs between AmQ-m vs AcQ-m and AmQ-c vs
   AcQ-c. d The numbers of DEGs and the up-down-regulation relationship in
   the same two comparisons using the A. cerana reference genome. e The
   DEGs in AmV vs AcV that showed the same expression trend using either
   A. mellifera or A. cerana as a reference genome. f The DEGs in AmQ vs
   AcQ that showed the same expression trend using either A. mellifera or
   A. cerana as a reference genome.

   AmV vs AcV common expression trend DEGs were enriched in Fatty acid
   degradation, mTOR, Glycerolipid metabolism and other pathways
   (Supplementary Fig. [85]2a). AmQ vs AcQ common expression trend DEGs
   were enriched in Wnt, mTOR, Noth, Glycerophospholipid metabolic
   pathway, and others (Supplementary Fig. [86]2b).

Differential open chromatin regions in queen ovaries between the two bee
species

   To further reveal the distribution of peaks in the genomes, 78,740 and
   42,324 peaks in AmV-m (A. mellifera virgin queens using the A.
   mellifera reference genome) and AmQ-m (A. mellifera laying queens using
   the A. mellifera reference genome) respectively were annotated to the
   genome of A. mellifera (Fig. [87]3a, [88]b), and there were 22,426 and
   17,826 differential peaks in the comparisons AmV-m vs AcV-m and AmQ-m
   vs AcQ-m (Fig. [89]3c). While 95,834 and 55,140 peaks in AcV-c (A.
   cerana virgin queens using the A. cerana reference genome) and AcQ-c
   (A. cerana laying queens using the A. cerana reference genome)
   respectively were annotated to the genome of A. cerana (Fig. [90]3d,
   e). The peaks were mainly distributed in the promoter region and intron
   region, which accounted for more than 54% of the total peaks. There
   were 19,288 and 15,383 differential peaks in the comparisons AmV-c vs
   AcV-c and AmQ-c vs AcQ-c (Fig. [91]3f). There were more up-regulated
   peaks in A. mellifera (AmV, AmQ) under both reference genomes
   (Fig. [92]3c, f). The number of differential peaks distributed in
   intron and promoter regions were predominant (Fig. [93]3g, h).

Fig. 3. ATAC-seq data analysis results.

   [94]Fig. 3
   [95]Open in a new tab

   a, b The distribution of peaks in AmV and AmQ in the genome of A.
   mellifera. c The numbers of differential peaks and the
   up-down-regulation relationship in the two comparisons (AmV vs AcV, AmQ
   vs AcQ) using the A. mellifera reference genome. d, e The distribution
   of peaks in AcV and AcQ in the genome of A. cerana. f The number of
   differential peaks and the up-down-regulation relationship in the same
   two comparisons using the A. cerana reference genome. g The
   distribution of differential peaks in the three comparisons AmV-m vs
   AcV-m, and AmQ-m vs AcQ-m. h The distribution of differential peaks in
   the three comparisons AmV-c vs AcV-c and AmQ-c vs AcQ-c.

Combined analysis of ATAC-Seq and RNA-Seq

   When comparing the differential peak-related genes with the DEGs, in
   AmV vs AcV, there were 6136 overlapped genes between ATAC-Seq and
   RNA-Seq, when using the A. mellifera reference genome. While there were
   6604 overlapped genes between ATAC-Seq and RNA-Seq, when using the A.
   cerana reference genome. Between these two sets of DEGs, 4411 genes
   overlapped (Fig. [96]4a). Analyzing these overlapped genes with the
   consistent trend genes in Fig. [97]2e, it was found that 861
   differential genes were shared (Fig. [98]4c). Of them, 88 genes were
   up-regulated and 84 genes were down-regulated in AmV in both ATAC-seq
   and RNA-seq under the A. mellifera reference genome (Fig. [99]4e), and
   the numbers were 131 and 85 under the A. cerana reference genome
   (Fig. [100]4g). Under the two reference genomes, there were 40
   consistently up-regulated genes and 16 consistently down-regulated
   genes in AmV in both ATAC-seq and RNA-seq (Fig. [101]4i). AmV vs AcV
   common trend DEGs were enriched in the fatty acid metabolism pathway
   and other pathways (Fig. [102]4j).

Fig. 4. Analysis the result of ATAC-Seq and RNA-Seq.

   [103]Fig. 4
   [104]Open in a new tab

   a, b Venn diagrams of DEGs and differential peak related genes in
   comparisons AmV vs AcV and AmQ vs AcQ using either A. mellifera or A.
   cerana as a reference genome. c, d Analyzing the common genes (a, b)
   with the consistent trend genes in Figs. [105]2e, f. e, f Using the A.
   mellifera reference genome, the RNA-seq and the ATAC-seq common DEGs of
   all comparison combinations in c, and d were statistically displayed.
   g, h Using the A. cerana reference genome, the RNA-seq and the ATAC-seq
   common DEGs of all comparison combinations in c and d were
   statistically displayed. i The panels e–h differential genes with the
   same trend in either reference genomes (the reference genome of A.
   mellifera or A. cerana). j, k The Kegg pathways of common trend DEGs.

   In AmQ vs AcQ, there were 5016 and 5599 overlapped genes between
   ATAC-Seq and RNA-Seq, when using the A. mellifera and A. cerana
   reference genomes, respectively. Of these, 3263 DEGs overlapped
   (Fig. [106]4b). Analyzing these overlapped genes with the consistent
   trend genes in Fig. [107]2f, it was found that 1037 DEGs were shared
   (Fig. [108]4d). Of them, 100 genes were up-regulated and 172 genes were
   down-regulated under A. mellifera reference genome (Fig. [109]4f), and
   the numbers were 131 up-regulated genes and 127 down-regulated genes in
   AmQ in both ATAC-seq and RNA-seq under the A. cerana reference genome
   (Fig. [110]4h). Under the two reference genomes, there were 48
   consistently up-regulated and 72 consistently down-regulated genes in
   AmQ in both ATAC-seq and RNA-seq (Fig. [111]4i). AmQ vs AcQ common
   trend DEGs were enriched on Wnt, mTOR and Glycerophospholipid metabolic
   pathway and other pathways (Fig. [112]4k).

Agxt2l, spr and nedd8 were essential for ovary development

   To identify candidate genes that led to differences in queen egg-laying
   between A. mellifera and A. cerana, from the 48 consistently
   up-regulated and 72 consistently down-regulated genes in AmQ vs AcQ in
   both ATAC-seq and RNA-seq using the dual species genomes, we chose 49
   genes that contained a differential peak in their promoter region
   (Fig. [113]5a) as candidates, and found three genes, agxt2l
   (LOC408817), spr (LOC725051) and nedd8 (LOC550745), were essential for
   ovary development. After RNAi of these three genes, their expression
   levels (Mann-Whitney U, P < 0.05), the number of ovarioles
   (Mann-Whitney U, P < 0.05) and body weight in the RNAi group
   (Mann–Whitney U, P < 0.05 and T-test, P < 0.05) (Fig. [114]5i–s) were
   all significantly reduced compared with the control group, except that
   the body weight showed no significant changes after RNAi of nedd8
   (T-test, t = −1.78, df = 38, P > 0.05) (Fig. [115]5t) (Supplementary
   Data [116]3).

Fig. 5. Screening and validation of candidate genes.

   [117]Fig. 5
   [118]Open in a new tab

   a The 49 genes that contain a differential peak in their promoter
   region. b–e Bees ovaries in the control group and RNAi agxt2l, spr and
   nedd8. f–h agxt2l, spr and nedd8 differed between AmQ and AcQ in all
   two omics assays. i–k Expression of agxt2l, spr and nedd8 in queen
   ovaries of both species. l–n Expression of agxt2l, spr and nedd8 in
   ovaries in RNAi injected queens and control. o–t Body weight and number
   of ovarioles in A. mellifera injected with siRNA of agxt2l, spr, nedd8
   and the control. u The glycerophospholipid metabolic pathway. Genes
   marked in red represent genes that have differences in gene expression
   be in AmQ vs AcQ. (we determined one side ovary because the other side
   was used for qRT-PCR. But we presented the total number here by timing
   it by 2 to be consistent with Fig. [119]1f).

Agxt2l was localized in nutrient cells

   Spatial orientation indicated that A. mellifera queen ovaries had more
   cells (nutrient, germ) than that of A. cerana queens (Fig. [120]6a).
   agxt2l was localized in nutrient cells (Fig. [121]6b) and expressed
   more in A. mellifera queen ovaries than in that of the A. cerana queens
   (Fig. [122]6c), details in supplementary Fig. [123]3. The results of
   DNA template preparation are in Supplementary Fig. [124]4.

Fig. 6. Spatial orientation of agxt2l expression in ovary.

   [125]Fig. 6
   [126]Open in a new tab

   a Cell type and proportion in ovary of A. mellifera queen and A. cerana
   queen. b Spatial orientation of agxt2l expression in ovary. c Results
   of in situ hybridization of agxt2l.

Brc-z1 regulates the transcription of agxt2l gene

   To explore the molecular mechanisms leading to the differential
   expression of the agxt2l gene between A. mellifera and A. cerana, we
   analyzed the genomic sequences corresponding to the differential peaks
   located in the promoter region of agxt2l in A. mellifera and A. cerana,
   and predicted two transcription factors whose binding sites differed
   between A. mellifera and A. cerana (Supplementary Fig. [127]5). To
   verify these three candidate transcription factors regulating the
   transcription of agxt2l, RNAi was performed in both honey bee species.
   It was found that after brc-z1 (LOC552255) RNAi, expression level of
   brc-z1 (N = 21, Z = − 4.49, P < 0.05) (Fig. [128]7a) and agxt2l
   (N = 21, Z = − 2.10, P < 0.05) (Fig. [129]7b), ovariole number
   (Mann–Whitney U, Z = −4.98, P < 0.05) (Fig. [130]7c) and body weight
   (T-test, t = −2.69, df  = 38, P < 0.05) (Fig. [131]7d) of A. mellifera
   queens in the RNAi group, were significantly lower than those in the
   control group (Supplementary Data [132]4). While the other genes and
   hsf (list the genes here) had no effect on the transcription of agxt2l.

Fig. 7. Motif analysis, transcription factor prediction and change of
phospholipid content after RNAi of agxt2l and brc-z1.

   [133]Fig. 7
   [134]Open in a new tab

   a, b Expression of brc-z1 and agxt2l injected with siRNA of brc-z1. c,
   d Body weight and number of ovarioles in A. mellifera injected with
   siRNA of brc-z1. e Change of phospholipid content after RNAi of agxt2l
   and brc-z1.

Agxt2l regulates phospholipid content in ovaries

   Phospholipid content in larval ovaries of agxt2l and brc-z1 RNAi groups
   were significantly lower than that of the control group
   (Kruskal–Wallis, H = 48.45, P < 0.001, P[adj] < 0.05), and the
   phospholipid content in agxt2l RNAi group was significantly higher than
   that in the brc-z1 RNAi group (Kruskal–Wallis, H = 48.45, P < 0.001,
   P[adj] < 0.05) (Fig. [135]7e) (Supplementary Data [136]5).

Agxt2l regulates egg-laying performance of queens

   The difference in egg-laying number before and after treatment was
   significantly greater in agxt2l RNAi group compared to the control
   group (T-test, t = 3.24, df = 16, P < 0.05) (Fig. [137]8)
   (Supplementary Data [138]7).

Fig. 8. Agxt2l regulates egg-laying performance of queens.

   Fig. 8
   [139]Open in a new tab

   The X-axis shows the experimental (injected agxt2l-siRNA solution) and
   control (injected negative control siRNA solution) groups, while the
   Y-axis shows the difference obtained by subtracting the
   pre-interference egg-laying quantity from the post-interference
   quantity.

Discussion

   The principal findings of this study are that we revealed the
   mechanisms underlying the higher egg-laying capability of A. mellifera
   over A. cerana. We identified agxt2l as a potential gene: its higher
   expression in trophocyte cells induces increased glycerophospholipid
   metabolism, accounting for the higher queen egg-laying rate in Apis
   mellifera compared to A. cerana. And agxt2l was in turn regulated by
   brc-z1 (LOC552255) (Figs. [140]5u,  [141]7,  [142]8,  [143]9, and
   Supplementary Fig. [144]5). In addition, brc-z1 promotes follicular
   opening. We determined that most of the enrichment pathways in
   egg-laying queens were related to oviposition or oocyte formation.
   Finally, we also verified that the egg-laying rate of A. mellifera
   queen is significantly higher than that of A. cerana queen, under the
   same environmental conditions (colony strength, season, and nectar
   resources).

Fig. 9.

   [145]Fig. 9
   [146]Open in a new tab

   Increased glycerophospholipid metabolism due to higher expression of
   agxt2l in trophocyte cell explains the higher queen egg laying rate in
   Apis mellifera than A. cerana.

   The egg-laying capability of A. mellifera queen and A. cerana queen
   have been studied many times^[147]33–[148]36 but no study controlled
   variables which could cause differences besides inherent capability. To
   rule out contribution due to other factors (such as
   seasons^[149]33,[150]35, population sizes^[151]34, regions^[152]36,
   etc.), we tested the egg-laying performance of the two species under
   the same conditions. We found that the number of eggs laid per
   milligram of body weight (Fig. [153]1f) and eggs per ovarioles
   (Fig. [154]1g) in A. mellifera queen was significantly higher than A.
   cerana queen, we hypothesize that the higher egg-laying rate in A.
   mellifera queen is mainly due to higher efficiency of each ovariole or
   higher maturation rate of oocytes in each ovariole.

   The differences in egg-laying rate might be due to differences in
   glycerophospholipid metabolic rates. In AmV vs AcV, the common trend
   DEGs were enriched in the fatty acid metabolism pathway (Fig. [155]5j).
   This pathway has been shown to play an important role in the
   development and reproduction of insects^[156]37–[157]41. In AmQ vs AcQ,
   the common trend DEGs were enriched in Wnt, mTOR signaling pathway and
   Glycerophospholipid metabolic pathway (Fig. [158]5k). Wnt signaling
   pathway has been found to be involved in the development of
   reproductive system, such as follicular development, ovulation, and
   luteinization^[159]42. mTOR and glycerolipid metabolism pathways
   indirectly promote each other. mTOR signaling has been shown to promote
   reproduction and lifespan in C. elegans^[160]43,[161]44, female D.
   melanogaster^[162]45 and female mice^[163]46. The first step of the
   Glycerophospholipid metabolic pathway is mediated by mTOR signal
   transduction^[164]47, which suggests that Glycerophospholipid metabolic
   pathway may be related to egg-laying difference between the two
   species. In addition, in honey bees, approximately 30–40% of the dry
   weight of an egg is consisted of lipids^[165]48. Glycerophospholipids
   are important lipid components in the queen ovaries and play crucial
   roles in the membrane structure and cell functions^[166]49, which
   suggests that glycerophospholipid metabolic rate is the reason for the
   difference in egg-laying rate.

   Agxt2l, spr and nedd8 promotes ovary development. We found that agxt21,
   spr and nedd8 were expressed higher in A. mellifera queen than in A.
   cerana, and they enhanced queen ovary development (Fig. [167]5c–e).
   Previous studies indicated that spr enhances egg-laying and female
   fecundity in Bactrocera dorsalis (Diptera: Tephritidae)^[168]50, D.
   melanogaster^[169]51 and Spodoptera litura^[170]52. Zebrafish nedd8
   facilitates ovarian development and the maintenance of female secondary
   sexual characteristics via suppression of androgen receptor
   activity^[171]53. This suggests that in honey bees, agxt2l, spr and
   nedd8 gene may be the candidate genes involved in differential
   egg-laying capacity between the two sister species.

   The rate of glycerophospholipid metabolic is promoted by agxt2l, which
   itself is enhanced by brc-z1. agxt2l is located in glycerophospholipid
   metabolic pathway and belongs to a subfamily of proteins called “Class
   II transaminases”, and is a modulator of phospholipid
   metabolism^[172]54. Liang et al. found that phospholipids promote
   ovarian development in shrimps^[173]55. In addition, in pregnant
   rabbits^[174]56, pigs^[175]57 and sheep^[176]58, phospholipids in the
   ovaries reach high levels during ovulation and the levels are
   maintained throughout pregnancy. These results indicate that
   phospholipids are related not only to ovarian development but also to
   embryo formation. We found that brc-z1 bound to the agxt2l promoter
   region (in both A. mellifera and A. cerana.) and positively regulated
   agxt2l. The broad-complex (br-c) is an early ecdysone hormone response
   gene that encodes four zinc finger protein subtypes (Z1, Z2, Z3, and
   Z4). In Blattella germanica, these four zinc finger isomers all play an
   important role in embryogenesis^[177]59. When BR-C is mis-expressed,
   the follicular cells in the front of the oocyte do not migrate forward
   properly^[178]60. In D. melanogaster, brc-z1 is involved in the
   formation of egg and ovum^[179]61–[180]63 and is specifically expressed
   in nurse cells, somatic follicle cells^[181]64, and egg chamber^[182]59
   which indicate that brc-z1 not only regulates agxt2l but also affects
   ovary development.

   Increased glycerophospholipid metabolism due to higher expression of
   agxt2l in trophocyte cell explains the higher egg-laying rate in A.
   mellifera queens. We found that agxt2l is highly expressed in
   trophoblast cells (Fig. [183]6b) and the proportion of trophoblast
   cells in A. mellifera queens is higher than that in A. cerana. In bees,
   trophoblast cells contain a large number of coarse endoplasmic
   reticulum and Golgi cytoplasm, mitochondria, lipids, carbohydrates and
   proteins^[184]65–[185]68. The presence of these organelles indicates
   that trophoblast cells are capable of synthesizing, absorbing,
   metabolizing, and storing substances, especially those derived from
   hemolymph^[186]69. After interfering with agxt2l, glycerophospholipid
   metabolism of trophoblast cells was limited, and lipid synthesis/
   transport was decreased (Fig. [187]9), thus perhaps resulting a delay
   (or reduction) in egg formation. Based on the present results, we
   provide evidence that the rate of lipid synthesis controlled by agxt2l
   is the main reason why the A. mellifera queens have higher egg-laying
   capability than the A. cerana queens.

   Furthermore, evolutionary pressure disparities likely underlie the
   variation in reproductive phenotypes. Apis mellifera ligustica is
   native to the Apennine Peninsula in the central Mediterranean, and Apis
   cerana cerana is mainly distributed in the subtropical humid area of
   the southern Longnan River Valley, the warm temperate humid area of the
   northern Longnan, the temperate semi-humid area of the southern
   Longzhong, and the temperate semi-arid area of the northern Longzhong,
   east of the Wushaoling Mountains in China. The monsoon climate in
   China, compared to Mediterranean climates, is characterized by high
   diversity, more acute changes, and a more scattered nectar flow. There
   are also large differences among different regions in China. All of
   these factors favor the phenotype for more flexible reproduction to
   achieve better resource management and faster responses to changing
   environments (including swarming and absconding), which increases the
   chance of survival for the colony under a changeable climate.

Materials and methods

Queen egg-laying capacity

   A. mellifera ligustica colonies (Jiangxi Agricultural University,
   28.46°N, 115.49°E, Nanchang, China) and A. cerana cerana colonies
   (Yifeng County, 28.26°N, 114.39°E, Jiangxi, China) (N = 10 each
   species) with newly laying queens (approximately 2 months old) were
   randomly selected for this experiment. The distance between the two
   places is 144 km. The queen of each colony was restricted to laying
   eggs on a comb for 24 h. The number of eggs on the comb was counted by
   taking digital photos (Canon EOD 70D) and using Image J (v1.54 d). This
   egg-laying experiment was repeated twice. We then fixed the queens on
   wax dish (N = 10 for each species) and dissected the ovaries. Ovary
   sections were prepared according to the method reported by Yi et
   al.^[188]70. Then, the number of ovarioles from both ovaries was
   counted under a dissection microscope (40x Nikon DS. Fi3). The
   differences in egg number and ovariole number between the two species
   were analyzed using Mann–Whitney Test and T-test in SPSS Statistics
   version 25.0.

Insects and sample collection

   A. mellifera ligustica colonies and A. cerana cerana colonies (N = 14
   colonies each species) in Jiangxi Agricultural University (28.46°N,
   115.49°E), Nanchang, China, were used in the rest of this study. For A.
   mellifera, the queen of each colony was confined on a comb to lay eggs
   for six hours (8:00 am-2:00 pm), and then the queen was separated from
   the experimental comb using a queen excluder. After egg hatching, the
   larvae were transferred to queen cells. On the tenth day after
   hatching, twelve queenless colonies were prepared as nucleus colonies.
   In the afternoon of the 11th day, some queen cells were placed in the
   queenless colonies to obtain laying queens, and some were placed in
   queen cages to obtain virgin queen samples. Seven ovary samples of the
   newly emerged virgin queen were collected. After the queens had mated,
   seven ovary samples of the newly laying queens (approximately two
   months old) were collected. The sampling methods of A. cerana were the
   same as A. mellifera.

   A total of 28 pairs of ovaries were collected from four types of
   queens, these were virgin queens and newly laying queens of both bee
   species (N = 7 for each type).

RNA extraction, transcriptome sequencing, and data analysis

   Total RNAs were extracted from the ovary samples using Trizol Reagent
   (Ambion/Invitrogen, USA). RNA integrity and total amount were detected
   by Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). All RNA
   samples that passed quality control were used to construct the
   library^[189]71. The library preparations were sequenced on an Illumina
   Novaseq platform (performed by Novogene Co., Ltd. (Beijing, China) and
   150 bp paired-end reads were generated. From the sequencing data, reads
   containing just adapters and undetermined bases as well as low-quality
   reads (Qphred≤20 reads with more than 50% of the total read length)
   were filtered out.

   The cleaned reads were aligned to the A. mellifera genome version
   Amel_HAv3.1(ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/254/395/GCF_003254
   395.2_Amel_HAv3.1/) or to the A. cerana reference genome
   ([190]ftp://ftp.ncbi.nlm.nih.gov/ assembly/ GCA_020141525.1) using
   HISAT2 software^[191]72. The expression level was calculated by
   Fragments Per Kilobase of transcript per Million (FPKM) base pairs,
   with standardized annotation of the length of the transcript and the
   total number of reads aligned with the transcriptome. Analysis of
   differential expression was performed using the DESeq2 R package
   (1.32.0). The p-value was corrected for multiple comparisons with a
   false discovery rate (FDR) < 0.05^[192]73–[193]75. Genes with padj≤0.05
   were defined as differentially expressed genes.

   After that, all DEGs were mapped to the Gene ontology (GO) database for
   enrichment analysis of functional significance. Finally, pathway
   enrichment analysis was performed in the Kyoto Encyclopedia of Genes
   and Genomes (KEGG) database^[194]76.

Nuclear extraction, ATAC-seq, and data analyses

   The ATAC-seq was performed according to Buenrostro et
   al.^[195]77–[196]79 and performed by Novogene Co., Ltd. (Beijing,
   China). Cell nuclei were extracted from the ovary samples and the
   nuclei particles were re-suspended in the Tn5 transpotase reaction
   mixture. The transposition reaction was incubated at 37 °C for 30 min.
   After transposition of Novogene, equimolar adapter1 and adapter2 were
   added to amplify the library by PCR. Then, the PCR products were
   purified by AMPure beads, and the quality of the libraries was
   evaluated by Qubit. The libraries were sequenced on an Illumina Novaseq
   platform.

   Index of the reference genome was built using BWA (v0.7.12) and clean
   reads were aligned to the reference genome using BWA mem, by default
   parameters. Peaks with q-value threshold of 0.05 were used for all data
   sets, peak-related genes were identified and annotate by
   ChIPseeker^[197]80. The peaks with fold change more than
   2(|log[2]fc | >1) were considered as differential peaks.

Genome homology alignment and data analyses

   The homologs between A. mellifera and A. cerana were identified using
   Blast software (2.2.26) with e value < 1e − 10. In the process of
   RNA-seq and ATAC-seq analyses, the homologous genes with the same trend
   were selected, and the differential genes and differential peaks
   showing difference when using the two sets of genomes as references