Extent and complexity of RNA processing in honey bee queen and worker caste development.

The distinct honeybee (Apis mellifera) worker and queen castes have become a model for the study of genomic mechanisms of phenotypic plasticity. Here we performed a nanopore-based direct RNA sequencing with exceptionally long reads to compare the mRNA transcripts between queen and workers at three points during their larval development. We found thousands of significantly differentially expressed transcript isoforms (DEIs) between queen and worker larvae. These DEIs were formatted by a flexible splicing system. We showed that poly(A) tails participated in this caste differentiation by negatively regulating the expression of DEIs. Hundreds of isoforms uniquely expressed in either queens or workers during their larval development, and isoforms were expressed at different points in queen and worker larval development demonstrating a dynamic relationship between isoform expression and developmental mechanisms. These findings show the full complexity of RNA processing and transcript expression in honey bee phenotypic plasticity.

Introduction

Phenotypic plasticity is a phenomenon where the same genotype showed distinct phenotypes if raised in different conditions, which is the most interesting aspect of development (Agrawal, 2001; Borges, 2005). Phenotypic plasticity is particularly dramatic in eusocial insects, and the honeybee (Apis mellifera) has become a valuable system for exploring this phenomenon. Honeybee queens and workers both develop from fertilized (diploid) eggs, but differences in nutrition in their larval diets determines their diverging developmental trajectories that result in different phenotypes (Asencot and Lensky, 1984; Brouwers, 1984; Maleszka, 2008; Spannhoff et al., 2011; Mao et al., 2015). The queen is the colony’s only fecund female with well-developed ovaries producing around 1500 eggs per day, whereas workers are facultatively sterile females that perform all the essential tasks in the colony such as brood care, foraging, defense, construction, and cleaning (Winston, 1991).

Prior foundational studies have shown how the process of phenotypic plasticity is epigenetically determined (Chen et al., 2012; Foret et al., 2012; Cameron et al., 2013; Guo et al., 2016; He et al., 2019; Wojciechowski et al., 2018). Studies have compared the transcriptomes between honeybee queens and workers using short-length RNA sequencing (RNA-Seq) (Chen et al., 2012; Cameron et al., 2013; He et al., 2019) or cDNA microarrays (Yamazaki et al., 2006; Barchuk et al., 2007). Thousands of genes were differentially expressed between queens and workers during their development (Chen et al., 2012; Cameron et al., 2013; He et al., 2019). Consistently a few key signaling pathways have been implicated as important in queen/worker differentiation. These include target of rapamycin (TOR), the fork head box O (FoxO), Notch, Wnt, insulin/insulin-like signaling (IIS) pathways, mitogen-activated protein kinase (MAKP), Hippo, and transforming growth factor beta (TGF-beta) (Patel et al., 2007; Mutti et al., 2011a; Duncan et al., 2016; Xiao et al., 2017; Yin et al., 2018; Wang et al., 2021). In addition, key genes such as Juvenile hormone esterase precursor (Jhe), ecdysone receptor (Ecr), Vitellogenin (Vg), and Hexamerin 70b (Hex70b) major royal jelly proteins (Mrjps) also participate in the determination of queen-worker fate (Barchuk et al., 2007; Buttstedt et al., 2013; Mello et al., 2014). Moreover, epigeneic modifications such as DNA methylation, microRNA, and alternative chromatin states have been shown to be involved in the regulation of honeybee queen-worker dimorphism (Foret et al., 2012; Ashby et al., 2016; Guo et al., 2016; Wojciechowski et al., 2018). The latest evidence showed that m6A methylation on RNA also participates in the genomic regulation leading to queen-worker dimorphism (Bataglia et al., 2021; Wang et al., 2021).

These studies have gone a long way to establish the honeybee as a model for genomic analyses of phenotypic plasticity, but thus far, we know rather little about how RNA processing might be involved in the process. RNA processing determines the mRNA coding sequence and can lead to different isoforms of a gene. Distinct isoforms can perform different biological functions, and these can be important for animal phenotypic plasticity (Marden, 2008). Alternative transcripts produced from the same gene can differ in the position of the start site, the site of cleavage and polyadenylation, and the combination of exons spliced into the mature mRNA (Parker et al., 2020). Poly(A) tails are also very important for RNA processing. Length of the Poly(A) tail is correlated with translational efficiency and regulates transcript stability (Workman et al., 2019; Roach et al., 2020).

Alternative splicing (AS) of genes diversifies the transcriptome and increases protein coding capacity through the production of multiple distinct isoforms (Modrek and Lee, 2001; Reddy et al., 2013). AS is involved in many biological processes in plants and animals and particularly responses and adaptation to changing environments (Modrek and Lee, 2001; Reddy et al., 2013). Alternative splicing of pre-mRNA has been associated with phenotypic plasticity in insects such as the bumble bee Bombus terrestris and the pea aphid Acyrthosiphon pisum (Grantham and Brisson, 2018; Price et al., 2018). Different alternative transcripts are also involved in honeybee caste differentiation, and alternative splicing plays an important role in queen-worker differentiation (Aamodt, 2008; Mutti et al., 2011b). These studies have focused on the effects of AS on the regulation of gene expression. Thus far, none have investigated how AS determines isoform expression genome wide in honeybees.

So far, honeybee RNA-Seq studies have been based on short-length RNA sequencing, which is limited to detect full-length isoforms or alternative transcripts of a gene. The new direct RNA sequencing technology (DRS) based on the Oxford Nanopore technology (ONT) directly sequences the native full-length mRNA molecules, which offers detailed information on the mRNA and RNA modifications (Garalde et al., 2018; Harel et al., 2019; Workman et al., 2019; Parker et al., 2020; Roach et al., 2020; Zhang et al., 2020). It avoids biases from reverse transcription or amplification and yields full-length, strand-special mRNA (Garalde et al., 2018; Harel et al., 2019; Zhang et al., 2020). This method allows a genome-wide investigation of RNA processing [e.g., alternative splicing and poly(A) tails].

Our objective in this study was to deploy this DRS technology to examine differences in RNA processing between the worker and queen honeybee castes. This allows us to investigate fundamental questions on caste-specific transcriptome patterns in honeybees: the variety of full-length isoforms and how they expressed differentially in queens and workers and the extent and complexity of the RNA processing in honeybee development.

We sampled honeybee queens and workers at different points in their larval development to assess the process of developmental canalization into the different phenotypes. Our results reveal extremely prevalent and highly complex differences in RNA processing between workers and queens. These differences changed across larval development showing a complex dynamical interaction between the epigenomic programs and the developmental pathways

Results

Quality of direct RNA sequencing data

Twelve libraries were generated and sequenced by ONT direct RNA sequencing. The number of clean reads of each sample was 5,416,798 ± 1,510,975 (Table 1). The N50 read length was 1,297 ± 30 nucleotides (nt), and the read quality was 10.63 ± 0.13 (Table 1), indicating a high integrity of the sequenced RNAs. 98.13 ± 0.79% clean reads were successfully mapped to the A.pis mellifera reference genome (Amel HAv 3.1) (Table 1). Moreover, the total number of genes detected from DRS and RNA-seq was equal (9,592 vs 9838), and over 88.8% of them were overlapped (Figure S1). We also mapped the transcript isoforms to honey bee reference genes and showed a great quality of percentile of the gene body (5′-3′) in each DRS sample (Figure S2). These results again confirmed a high and reliable integrity of the sequenced RNAs. The Pearson correlation coefficients between two biological replicates of each group were all over 0.85 (Figure S3) and the sample cluster tree also showed that two biological replicates of each group were clustered together (Figure S4), reflecting a good repeatability of biological replicates.

Table 1

Overview of direct RNA sequencing

Samples	2Q	2W	4Q	4W	6Q	6W	Average
Clean read numbers	6,531,006 ± 1,043,238	6,103,488 ± 14,688	6,506,243 ± 620,747	5,007,500 ± 1,696,509	3,548,493 ± 2,108,736	4,804,061 ± 1,564,551	5,416,798 ± 1,510,975
Average length (nt)	1,030 ± 46	1,012 ± 15	1,003 ± 7	998 ± 42	1,053 ± 23	1,002 ± 36	1,016 ± 31
N50 length (nt)	1,316 ± 57	1293 ± 11	1,286 ± 13	1,279 ± 52	1,311 ± 39	1,296 ± 14	1,297 ± 30
Max read length (nt)	20,044 ± 4,581	17,248 ± 987	24,924 ± 7,011	17,955 ± 2,212	17,548 ± 398	14,579 ± 1,298	18,716 ± 4,273
Average read quality	10.66 ± 0.07	10.73 ± 0.03	10.40 ± 0.19	10.64 ± 0.03	10.73 ± 0.06	10.65 ± 0.01	10.63 ± 0.13
Alignment % identity	98.58 ± 0.32	98.44 ± 0.03	98.43 ± 0.28	98.29 ± 0.13	97.63 ± 1.10	97.45 ± 1.78	98.13 ± 0.79
Uni mapped genes	8,736 ± 64	8,658 ± 0	8,544 ± 159	8,323 ± 125	8,204 ± 654	8,592 ± 74	8,509 ± 286
No. of isoforms	23,518 ± 764	22,896 ± 61	22,382 ± 899	20,144 ± 501	19,426 ± 4,156	20,800 ± 121	21,528 ± 2,040

The 2Q, 2W, 4Q, 4W, 6Q, and 6W are 2-day queen larvae, 2-day worker larvae, 4-day queen larvae, 4-day worker larvae, 6-day queen larvae and 6-day worker larvae, respectively. Data are represented as mean ± SD and each group has two biological replicates.

Differences of isoform and gene expression in queen-worker comparisons

Our results showed genome-wide differentially expressed isoforms in honeybee caste differentiation. Significantly differentially expressed transcript isoforms (DEIs) and genes (DEGs) between queen and worker larvae were identified. We identified 662, 1855 and 1042 DE-Is in 2d, 4d, and 6d queen-worker larval comparisons, respectively, which were notably more than the DEGs (281, 1369, and 645, respectively, see Figure 1A, detailed information see Tables S1, S2, S3, S4, S5, and S6). To compare the differences between DE-Is and DE-Gs, we firstly mapped DE-Is to reference genes (DEIGs) and showed that only 271 (34.72%), 977 (53.45%), and 410 (38.32%) genes overlapped between DEIGs and DE-Gs in 2d, 4d, and 6d queen-worker comparisons, respectively (Figures 1B–1D). This shows a considerable difference between DEIs and DEGs in queen-worker differentiation.

Figure 1

Comparison of isoform and gene expression, poly(A) length between queen and worker larvae

See also Tables S1, S2, S3, S4, S5, S6, S7, S8, S9, S13, S14, S15, S16, and S17.

(A) Numbers of differentially expressed isoforms (DEIs) and genes (DEGs), different poly(A)-length related isoforms between queens and workers at three larval stages (day 2, day 4, and day 6). The up bars are DEIs and DEGs length highly expressed in 2d, 4d, and 6d queen larvae as well as isoforms with longer poly(A) tails in queen larvae, and the down bars are that highly expressed or with longer poly(A) tails in worker larvae.

(B) The venn diagram of DEIGs (DEI mapped genes) and DEGs from 2d queen-worker comparison.

(C) The venn diagram of DEIGs and DEGs from 4d queen-worker comparison.

(D) The venn diagram of DEIGs and DEGs from 6d queen-worker comparison.

We identified DEIs (24, 54, and 48 in 2d, 4d, and 6d comparisons, respectively) that are involved in several key KEGG signaling pathways for honey bee caste differentiation, such as mTOR, Notch, FoxO, Wnt, IIS, Hippo, TGF-beta, and MAPK (see Tables S1, S2, and S3). The number of DEGs enriched in these signal pathways was 7, 34, and 33 in 2d, 4d, and 6d, respectively (Tables S4, S5, and S6). Other DEIs mapped to key genes such as Jhe, Mrjps, Ecr, Vg, and Hex70b (see Tables S1, S2, and S3). The number of isoforms with significantly different poly(A) length in 2d, 4d and 6d comparisons showed a negative trend as DEIs and DEGs, with more in queens (Figure 1A, detailed information see Tables S7, S8, and S9) and less in workers.

DEIs and DEGs were enriched in five key KEGG signaling pathways (mTOR, Notch, FoxO, Wnt, and IIS) which have been previously shown to be participating in queen-worker differentiation. DEIs were mapped to 35 proteins in these five pathways, whereas DEG mapped onto just 17 proteins (Figure 2A). In detail, 23 proteins were uniquely mapped by DEIs, whereas only 5 were uniquely mapped by DEGs. 12 proteins mapped to both DEI and DEG (Figure 2A). We also provided detailed information about the expression of DEIs and DEGs enriched in these five key pathways (Figure 2B). Clearly, more isoforms were significantly differentially expressed in these five key pathways compared to DE-Gs, and the differentially expressed isoforms of each DEIG are presented in Figure 2B.

Figure 2

Expression and enrichment of DEIs and DEGs in five KEGG pathways

See also Tables S1, S5, S6, S16, and S17.

(A) DEIs and DEGs enriched in key steps of five KEGG pathways. Five KEGG pathways were shown with transparent boxes and the names of these pathways were red marked. Green boxes inside of transparent boxes are DEI enriched proteins; yellow boxes are DEG enriched proteins; green-yellow mixed boxes are both DEI and DEG enriched proteins; blank boxes are non DEI and DEG enriched proteins.

(B) Expression of DEIs and DEGs in the five key pathways. The left are pathway and protein names, and the right are isoform expression and names. Expression of isoforms and genes are presented with their log2 TPM values and shown with color scales. DEIs or DEGs are marked with “∗” in the middle of the boxes.

These findings suggest that caste-specific isoform expression provided more accurate and detailed information of the real transcriptome differences between honeybee queen and worker castes than measurements of differences in gene expression.

Uniquely expressed isoforms between queens and workers

We identified 187, 357 and 364 uniquely expressed isoforms in queen or worker larvae at 2, 4aysd, 4d, and 6d age, with more occurring in queens than workers (Figure 3A). The number of uniquely expressed isoforms in queens reached a maximum at the 4d stage, whereas in workers it reached a maximum in the 6d sample. This could suggest that the queen developmental pathway diverges more quickly and at an early larval stage, whereas the developmental pathway of workers diverges slightly later.

Figure 3

Uniquely expressed isoforms in queen and worker larvae

See also Tables S10, S11, S12b, S16 and S17.

(A) Numbers of uniquely expressed isoforms in queen larvae (up bars) or worker larvae (down bars) at 2d, 4d and 6d stages. The red bars are numbers of uniquely expressed isoforms enriched in eight key KEGG signaling pathways (mTOR, IIS, FoxO, Notch, Wnt, MAPK, Hippo, and TGF-beta). The purple bars are numbers of uniquely expressed isoforms that are key genes reported in previous studies (details see Tables S10, S11, and S12).

(B) The heat map of expression of uniquely expressed isoforms which are shown in red bars and green bars in (A). The red color blocks in heat map are uniquely expressed isoforms in queen larvae at 2d, 4d and 6d age, whereas the dark green blocks are uniquely expressed isoforms in worker larvae. Other color scales are the log2 fold change values of isoform expression between queen and worker larvae at 2d, 4d, and 6d stages.

Uniquely expressed isoforms were enriched into eight key KEGG signaling pathways such as mTOR and IIS for queen-worker differentiation as well as some key genes such as DNA methyltransferase 3 (Dnmt3) and Ecr (Figure 3, details see Tables S10, S11, and S12). This indicated that there are isoforms which are likely to be involved in queen-worker differentiation that are uniquely expressed in queens or workers during their development.

DEIs and alternative splicing

We identified 21,574 isoforms in each sample on average, and these isoforms could map to averagely 8509 honeybee reference genes (Table 1). This suggested that the honey bee genome produces notably more isoforms than genes. DRS technology allows detecting the splicing events causing an isoform. We reconstructed the complex splicing events as well as the DEIs resulting from the alternative splicing. We then define a DEI from a single gene, which is differentially spliced compared to the longest isoform of this gene as a DS-DEI. The majority of DEIs (73.26%, 67.98, and 63.53% of 2d, 4d, and 6d DEIs) were DS-DEIs, and most of them (71.13%, 66.16%, and 65.56% in 2 d, 4d, and 6d) were generated by at least two types of splicing events in a single isoform. A combination of RI, A5, and A3 in a single DEI was the most common type (Figure 4). This indicates that the formation of transcripts in honeybee development is more complex than previously known.

Figure 4

The alternative splicing events (AS) in DEIs of three comparisons

See also Figure S6.

(A) DEIs of 2d comparison containing AS events. The colorful bars in the left bottom diagram are the different AS patterns and the number of AS events in DEIs. Seven different AS patterns are 3'splice site (A3), 5'splice site (A5), First exon (AF), Last exon (AL), Retained intron (RI), Skipping exon (SE), and Mutually exclusive exon (MX). The pink vertical bars in the top of the right diagram are the numbers of DEIs that were spliced by a particular combination of AS forms. The black spots in a column mean DEIs that were spliced by a combination of related AS forms. Same in (B and C).

(B) DEIs of 4d comparison containing AS events.

(C) DEIs of 6d comparison containing AS events.

Of these DEIs, 9.82% (65), 5.71% (106), and 9.50% (99) of DEIs from the 2d, 4d, and 6d comparisons were produced by significantly different AS events when comparing workers and queens (Figure S5, Tables S13, S14, and S15). Similarly, many of them (46.15%, 42.45, and 45.45% in 2d, 4d, and 6d comparisons, respectively) contained at least two types of significantly different AS events (Figure S6, Tables S13, S14, and S15). Consequently, a part of DEIs were significantly differentially spliced.

Correlation between poly(A) length and expression of DEIs and DEGs

We showed a negative correlation between isoform expression and their poly(A) tail length (Figure S7). Expression of DEIs (Figure 5) and DEGs (Figure S8) was also negatively correlated with their poly(A) lengths. It is believed that poly(A) tails are involved in the regulation of isoform expression. More interestingly, the poly(A) lengths were positively correlated with log2 fold change values of DEIs, which was stronger than that of DEGs (Figures 5D–5F). This firstly showed that poly(A) lengths were more closely related to the biased expression of DEIs than DEGs.

Figure 5

Correlation between poly(A) length and expression of DEIs and DEGs

See also Figures S7 and S8 and Tables S7, S8, and S9.

(A) Correlation between poly(A) lengths and the expression of DEIs of 2d queen-worker comparison. Expression of each DEI was the log10 TPM value. Same in (B and C).

(B) Correlation between poly(A) lengths and DEIs of 4d comparison.

(C) Correlation between poly(A) lengths and DEIs of 6d comparison.

(D) Correlation between poly(A) lengths and log2 fold change values of DEIs and DEGs of 2d comparison. Same in (E and F).

(E) Correlation between poly(A) lengths and log2 fold change values of DEIs and DEGs of 4d comparison.

(F) Correlation between poly(A) lengths and log2 fold change values of DEIs and DEGs of 6d comparison.

Correlation of isoform expression and poly(A) lengths for two key genes

Two key genes, Jhe and Ecr, which participate in queen-worker differentiation were selected as examples to show the correlation pattern of gene and isoform expression and its poly(A) length. Jhe had only one isoform and was negatively correlated with its poly(A) lengths (Figure 6A). The Ecr gene has 5 isoforms and we detected 4 isoforms (Figure 6B). Queen and worker larvae showed a more detailed difference in the expression of these 4 isoforms (see the middle of Figure 6B) than gene expression (see the left of Figure 6B). We selected the most differentially expressed isoform (Ecr.t4) to measure the correlation between its expression and poly(A) length. Similarly, the expression of Ecr.t4 isoform was also negatively correlated with their poly(A) lengths (Figure 6B).

Figure 6

The poly(A) lengths and isoform expression of two key genes

See also Tables S1, S2, S3, S16, and S17.

(A) Expression and related poly(A) lengths of Jhe isoform. The left is the exon coverage of clean reads (namely gene expression) of Jhe gene of 6 larval groups. The Jhe has only one isoform (Jhe.t1) and its structure is present (see its exons below marked with black color). The middle heat map is the isoform expression of Jhe.t1 (log2 TPM values) in 6 groups and presented with color scales. The right is the poly(A) length distribution of Jhe.t1 in 6 groups. The mean poly(A) length is presented in the middle of each violinplot.

(B) Expression and related poly(A) lengths of Ecr isoforms. The left is the exon coverage of clean reads (namely gene expression) of Ecr gene The Ecr gene has 4 isoforms and their expression is presented using heat map as above. The poly(A) length distribution of Ecr.t4 in 6 groups is presented on the right side; it was the most significantly differentially expressed isoform of the gene Ecr. The poly(A) of Ecr.t4 in 2dW was presented with a blank box, which means few clean data of poly(A) was detected because 2dW did not express Ecr.t4.

We selected two isoforms of the Ecr gene (Ecr.t1 and Ecr.t2) for TaqMan real-time PCR confirmatory experiment. The expression of these two isoforms in TaqMan real-time PCR (Figure 7) was mainly consistent with the DRS results, confirming a biased isoform expression in honey bee caste differentiation.

Figure 7

Verification of two isoforms of Ecr gene between 2d, 4d, and 6d queen and worker larvae by TaqMan real-time PCR

See also Table S18.

(A) The Ecr.t1 expression in three queen-worker comparisons (2 d, 4 d, and 6 d) from DRS. The expression of the isoform was the FPKM value. Data were presented as mean ± SE. Different letters in the top of bars indicate significant difference (p-values <0.05 and log2 Fold change values > 1.5 using DESeq2) and the same letter indicates no significant difference. Same in (B).

(B) The Ecr.t2 expression in three queen-worker comparisons (2d, 4d, and 6d) from DRS.

(C) The Ecr.t1 expression in three queen-worker comparisons (2d, 4d, and 6d) from TaqMan real-time PCR. Data were presented as mean ± SE. Different letters in the top of bars indicate significant difference (p < 0.05, one-way ANOVA test) and same letter indicate no significant difference. Same in (D).

(D) The Ecr.t2 expression in three queen-worker comparisons (2d, 4d, and 6d) from TaqMan real-time PCR.

Discussion

RNA processing acts as a major driver of animal phenotypic plasticity (Gingeras, 2007), but the complexity and characteristics of RNA processing underlying this phenotypic variability is not well understood. The present study has shown the extent of RNA processing in honeybee queen-worker differentiation. Differences in isoform expression revealed the true extent of the differences in the transcriptome between queens and workers during their development, providing more detailed and accurate information on mRNA expression compared to measures of gene expression. We also showed the flexibility of alternative splicing in queen-worker differentiation. Most DEIs were formed by more than one type of alternative splicing. Poly(A) length negatively correlated with the expression of DEIs demonstrating another part of the genomic mechanisms of honey bee caste differentiation.

Thousands of DEIs were identified between queens and workers, which were notably more than DEGs (Figure 1A). Many genes were not differentially expressed, but different isoforms from these genes were significantly differentially expressed (Figures 1B–1D, Tables S1, S2, and S3), demonstrating how DRS gives a very different perception of the transcriptome differences between developing worker and queen honey bees. DEIs are enriched into more proteins of five key KEGG signaling pathways such as mTOR and IIS than DEGs (Figure 2). The mTOR signaling pathway is the central component of a conserved eukaryotic signaling pathway, regulating cell and organismal growth in response to nutrient status and has been shown to be important for determining the queen-worker developmental paths (Colombani et al., 2003; Oldham and Hafen, 2003; Patel et al., 2007). The IIS and Notch signaling pathways have been shown to be involved in the development of honeybee ovaries (Mutti et al., 2011a; Duncan et al., 2016). More DEIs were identified in these key pathways than DEGs, indicating that DEIs could provide more nuanced differences of mRNA expression between honeybee queens and workers during their development. Therefore, compared to DE-Gs that conceal a more nuanced picture of transcription differences, DE-Is show the true nature of transcriptome differences between queens and workers. Consequently, DEIs identified in the present study revealed a more informative knowledge and precise difference in the transcriptome expression of honey bee caste differentiation compared to previous RNA-Seq results (Chen et al., 2012; Cameron et al., 2013; He et al., 2019; Yin et al., 2018).

Hundreds of uniquely expressed isoforms were identified in queens or workers. Many uniquely expressed isoforms were enriched in mTOR, IIS, FoxO, Notch, Wnt, MAPK, TGF-beta, and Hippo—the eight key KEGG-signaling pathways (Figure 3). Some uniquely expressed isoforms mapped to key genes such as DNA methyltransferase 3 (Dnmt3) and Ecr genes (Figure 3). Dnmt3 is a key driver of epigenetic global reprogramming controlling honeybee queen-worker fate (Kucharski et al., 2008). Here, we showed one isoform of Dnmt3 (Dnmt3.t3) was uniquely expressed in 4 days worker larvae, suggesting that this isoform of Dnmt3 may play an important role in the worker developmental pathway at this critical stage. For the Ecr gene, 2 days worker uniquely expressed Ecr.t3, whereas 2 days queen larvae uniquely expressed Ecr.t4. Different isoforms from the same gene were expressed in queens and workers, revealing the complexity of RNA processing in the honeybee queen-worker differentiation. We have documented many genes where isoform expression was not consistently biased toward either queens or workers throughout all stages of larval development. In some cases, an isoform was uniquely expressed in either queens or workers in one larval stage, but at other developmental stages the isoform occurred in both workers and queens (Figure 3B). Isoforms like LOC409856.t4 and LOC727227.t8 uniquely expressed in queens at 2d and 4d stages but oppositely expressed in workers at 6d stage (Figure 3B). This illustrates the dynamic relationship between RNA processing and the developmental process. A simple model might imagine different isoforms of a gene specific to either the queen or worker developmental pathway. Instead, we see isoforms featuring in both pathways but at different stages.

The number of uniquely expressed isoforms also showed a different trend in queens and workers. It reached a maximum in queens at the 4d larval stage but it peaked later in workers (Figure 3A). These different trends correlate with the known faster development of queens (Winston, 1991) and suggest the queen’s developmental pathway differentiates at an earlier larvae stage than the worker developmental pathway. Queen larvae at early stage receive significantly more food than worker larvae, and the nutrient contents in royal jelly also differ from worker jelly in terms of sugar, vitamins, proteins, acids, and microRNA (Asencot and Lensky, 1984; Brouwers, 1984; Maleszka, 2008; Mao et al., 2015). In Figure 3B, most isoforms enriched in mTOR signaling pathway which is a signaling pathway in response to nutrient status (Colombani et al., 2003; Oldham and Hafen, 2003) were uniquely expressed in queens at 2d stage and its number decreased at 4d and 6d stages, which supports our speculation.

We detected 21,528 ± 2,040 isoforms in each sample which could map to 8,509 ± 286 honey bee reference genes (Table 1). This is consistent with a previous study that 60.3% of honeybee multi-exon genes are alternatively spliced (Li et al., 2013). The majority of DEIs, contained AS events and over 65% of them had at least two types of AS patterns in a single isoform (Figure 4). Some DEIs contained more than five different AS types (Figure 4). A few DE-Is were differentially spliced between queens and workers by multiple AS patterns (Figure S3, Tables S13, S14, and S15). It means several distinct AS types together contributed to the construction of a single isoform (Figure 4), revealing a flexible system in honey bee caste differentiation.

Poly(A) tails are a regulator of translation and transcript stability (Lim et al., 2016; Tudek et al., 2018; Woo et al., 2018). Poly(A) tail length is not stable but dynamic and condition-dependent (Lim et al., 2016; Woo et al., 2018). This study showed that poly(A) tail lengths were strongly and negatively correlated with expression of isoforms in honeybees (Figure 5). The negative correlation between isoform expression and poly(A) length has been shown in the regulation of Caenorhabditis elegans development (Roach et al., 2020). Here we showed that poly(A) tails also participate in the regulation of isoform expression and therefore contributes to complexity of RNA processing in honey bee queen-worker differentiation. Moreover, we examined the relationship between poly(A) length and the degree of biased expression of isoforms or genes in queen-worker differentiation. Interestingly, the log2 fold change values of DEIs were positively correlated with the poly(A) lengths and their correlation were stronger than that of DEGs (Figures 5D–5F), suggesting that poly(A) tails play an important role in the RNA processing of honeybee caste differentiation. This also firstly showed that poly(A) tails had a stronger correlation with isoforms than genes in terms of insect phenotypic plasticity. Consequently, poly(A) tails possibly contribute to the determination of honeybee queen-worker dimorphism via regulating the isoform expression.

Environmental stimulus altering patterns of gene expression has been considered as a primary molecular mechanism of phenotypic plasticity (Schlichting and Simth, 2002; Aubin-Horth and Renn, 2009). However, one gene can produce different transcripts by RNA processing, and different transcripts have distinct biological functions that are vital for phenotypic plasticity (Marden, 2008). Therefore, increasing evidence revealed that the complex RNA processing may be a principal contributor to the phenotypic complexity (Gommans et al., 2008, 2010). The present study supports this hypothesis and showed a genome-wide different isoform expression in queen-worker differentiation, which was far more complex than the gene expression (Figures 1 and 2). Especially hundreds of isoforms were uniquely expressed in queens or workers (Figure 3). These uniquely expressed isoforms cannot be easily seen in gene expression but may play an important role in the determination of honey bee caste differentiation. We also showed a flexible alternative splicing generating a variety of isoforms (Figure 4 and S6) and the polyadenylation [poly(A) tails] negatively regulating isoform expression (Figures 5 and 6). All findings demonstrated a complexity of RNA processing in honey bee phenotypic plasticity.

In summary, RNA processing plays an important role in shaping honeybee phenotypic plasticity, but investigations of gene expression by short-length RNA-Seq fail to reveal the full complexities of RNA processing (Chen et al., 2012; Cameron et al., 2013; He et al., 2019; Yin et al., 2018). By using DRS, this study has shown the extent of differential isoform expression between workers and queens, flexibility of transcript splicing, and polyadenylation. It provides a more detailed understanding of the molecular mechanisms underlying the divergence of the queen and worker developmental paths. It also contributes to our understanding of the extent and complexity of RNA processing in animal phenotypic plasticity.

Limitations of the study

Although this study showed an extensive and complex RNA processing in honeybee caste differentiation including biased isoform expression, flexible alternative splicing, and polyadenylation, other RNA processing events like 5′ capping are limited to be detected by the DRS technology but may also be involved into the honeybee caste differentiation. In addition, to extract enough amount of total RNA for DRS, we used a whole body of a 4d and 6d queen or worker larva as a sample and collected twenty 2d larvae as a mixed sample because of its small body size. This amount of RNA used in DRS was much lower than that in RNA-Seq; however, we note that the isoform expression of honey bee queen and worker larvae might differ in their distinct tissues.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Prof. Zhi Jiang Zeng (bees1965@sina.com).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Insects

Three healthy honeybee colonies (Apis mellifera) each with a mated queen were used for this study. Each colony had ten frames with approximately 35,000 bees. These colonies were maintained at the Honeybee Research Institute, Jiangxi Agricultural University, Nanchang, China, according to the standard beekeeping techniques.

Method details

Larva sampling

The mother queen was caged onto a plastic worker-cell frame developed by Pan et al. (2013) to lay eggs into worker cells for 6 h. Afterwards, half of these eggs were removed into commercial plastic queen cells to rear queens before hatching. Queen and worker larvae at 48 h (2 days), 96hrs (4 days) and 144 h (6 days) after hatching were harvested and were flash-frozen in liquid nitrogen. In the DRS experiment, for the 2 days queen larvae (2Q) or 2 days worker larvae (2W), larvae were very small, therefore we collected twenty 2days queen larvae and mixed them into one sample. The 4d and 6d larvae were large enough, and one sample needed only one larva. Each group had two replicates from two colonies, therefore in total 12 samples were collected. Larval samples for short-length RNA-Seq were sampled the same as DRS, but with twice the number of larvae in each sample. We had three replicates for each of these 6 groups: 18 samples in total. Within a colony queens and workers developed from a same mother queen and eggs were laid at the same time. All sampled eggs were reared in a same colony. We compared queen and worker larvae sampled at the same time point: 2Q vs 2W, 4Qvs 4W, and 6Q vs 6W.

RNA preparation

Total RNA of 2d, 4d and 6d queen and worker larvae as extracted in accordance with the standard protocol of the of the Total RNA Kit I (R6834, Omega Bio-tek), respectively. Each sample contained 30 μg total RNA for DRS. RNA was cleaned and concentrated by NEBNext Poly(A) mRNA Magnetic Isolation Module (E7490S) according to the manufacturer’s instructions. From each sample we extracted 1 μg of total RNA for the RNA purity and concentration measurement using Nanodrop (Thermo Fisher Scientific). RNA integrity was checked using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA).The RNA samples with RNA Integrity Number (RIN) >8.5 and 2.0 < OD260/280 < 2.2 (Zhang et al., 2019) were used for DRS library preparation.

Direct RNA sequencing

Prepared RNA of each sample was used for a DRS library preparation using the Oxford Nanopore DRS protocol (SQK-RNA002, Oxford Nanopore Technologies). For reversed connector connection, 9 μL prepared RNA, 3 μL NEBNext Quick Ligation Reaction Buffer(NEB), 1 μL RT Adapter (RTA) (SQK-RNA002) and 2 μL T4 DNA Ligase (NEB) were mixed together and incubated under 25°C for 10 min. Afterwards, 8 μL 5× first-strand buffer (NEB), 2 μL 10 mM dNTPs (NEB), 9 μL Nuclease-free water, 4 μL 0.1 M DTT (Thermo Fisher) and 2 μL Super-Script III Reverse Transcriptase (Thermo Fisher Scientific, 18,080,044) were added into the above 15 μL reaction system and incubated under 50°C for 50 min then 70°C for 10 min. Reverse-transcribed mRNA was purified with 1.8∗Agencourt RNAClean XP beads and washed with 23 μL Nuclease-free water, and subsequently sequencing adapters were added using 8 μL NEBNext Quick Ligation Reaction Buffer, 6 μL RNA Adapter (RMX) and 3 μL T4 DNA Ligase. The mix was purified and washed again as above, and 75 μL RRB(SQK-RNA002)were added with 35 μL Nuclease-free water. This final reaction system was loaded into Nanopore R9.4 sequencing micro-array and sequenced for 48–72 h using PromethION sequencer (Oxford Nanopore Technologies). This was performed by Wuhan Benagen Tech Solutions Company Limited.

Preprocessing and alignments

To assess the read quality, raw reads were base called using Guppy software (version 3.2.6) (Fesenko et al., 2020) and adapter sequences were trimmed by Guppy as well. Low quality reads (Q-value < 7) and short-length reads (<50 nt) were filtered by Nanofilt (version 2.7.1) (Coster et al., 2018). The remaining clean reads were corrected using the Fclmr2 (version 0.1.2) with our short-length RNA-Seq data (Wang et al., 2018). Afterwards, clean reads from each library were aligned to the genome of honeybees Apis mellifera (Amel HAv3.1) using Minimap2 (version 2.17-r941) (Li, 2018). The alignment ratio of clean reads to the honeybee reference genes were calculated using Samtools (version 1.10) (Li et al., 2009) and showed in Table 1. After clearing redundant sequences using Flair collapse (version 1.5.0) (Tang et al., 2020) and Stringtie (version 2.1.4) (Pertea et al., 2015), the clean reads from each library were merged by Gffcompare (version 0.12.1) (Pertea and Pertea, 2020) to obtain the final isoforms. To test the RNA integrity of DRS samples, raw reads were used to calculate its distribution density on the gene body of honey bee genome by RSeQC software (4.0.0 version, Wang et al., 2012), and results were shown in Figure S2.

DEG and DEI calculation

Isoform expression was annotated to the honeybee reference genes and quantified using the selective-alignment-base model in Salmon software (version 1.4.0) (Patro et al., 2017). The raw sequence and annotation information of each isoform refer to Tables S16 and S17 respectively. Gene expression was measured according to standard methods as in He et al. (He et al., 2019) using DRS raw long reads rather than RNA-Seq data, and then Significantly differentially expressed genes (DEGs) between queen and worker larvae were identified as those with p values <0.05 and log2 Fold change values >1.5 using DESeq2 (version 1.26.0) (Love et al., 2014). Since DRS provides the expression of each isoform, therefore significantly differentially expressed isoforms (DEIs) between queen and worker larvae were identified using the same analysis (p values <0.05 and log2 Fold change values >1.5 using DESeq2). To compare the expression level of genes/isoforms, we computed TPM values (transcripts per kilobase of exon model per million mapped reads) based on reads counts as expression. Since it could not compare the DE-Gs and DE-Is directly, therefore, all DEIs were referenced to their related genes to obtain differentially expressed isoforms of genes (DEIGs). Sequences of the DE-Gs and DEIGs were blasted to the Swiss-Prot database, non-redundant protein (Nr) and non-redundant nucleotide sequence (Nt) database with a cut-off E-value of 10−5.

KEGG pathway enrichment analysis

Five KEGG pathways [TOR, FoxO, Notch, Wnt and IIS] which have been repeatedly associated with honeybee queen-worker differentiation (Patel et al., 2007; Mutti et al., 2011a; Duncan et al., 2016; Xiao et al., 2017; Yin et al., 2018; Wang et al., 2021) were selected and combined into a pathway net according to KEGG pathway database and a previous study (Wang et al., 2021). The statistical enrichment of DEIGs and DEGs in these KEGG pathways were analyzed using a hypergeometric test (Q-value < 0.05) using the KOBAS 2.0 software (Xie et al., 2011). Results were shown in Figure 3A. Expression of DE-Is and genes in these five pathways were also presented in Figure 3B.

Alternative splicing identification

The AS events of each comparison were analyzed using SUPPA2 (version 2.3) (Trincado et al., 2018). AS events were divided into seven patterns including 3’splice site (A3), 5’splice site (A5), First exon (AF), Last exon (AL), Retained intron (RI), Skipping exon (SE) and Mutually exclusive exon (MX). The PSI values of each AS events in DE-Is were calculated after filtering percent spliced in index (PSI) values with ‘not available’ (NA) or ‘0’ (Trincado et al., 2018). Significant differences in AS events were identified as its dPSI-value >0 and p value <0.05.

Poly(A) length estimation and correlation with expression of genes and isoforms

Poly(A) lengths of each read were calculated by the Nanopolish software (version 0.13.2) (Coster et al., 2018). The 3′ poly(A) tail lengths of each read were directly identified using the ionic current signals (Workman et al., 2019). Low-quality poly(A) tails was filtered by the “qc_tag” software in Nanopare sequencing platform, and only the PASS reads were used for further analysis. A previous study showed that poly(A) lengths were negatively correlated with isoform expression, therefore, the median of poly(A) length of each isoform was used for correlation and regression analysis with their TPM values (results were shown in Figure S5) in R package (version 4.0.3). We also evaluated the correlation between poly(A) tail lengths and DEIs or DEGs.

For this we used median of poly(A) length and TPM or log2Fold change values of DEIs or DEGs for correlation and regression analysis as above. Two key genes (Jhe and Ecr) for the determination of queen-worker differentiation were selected as an example to show the correlation between expression of isoforms and poly(A) lengths. Results were shown in Figure 6.

Short-length RNA sequencing

Total RNA of each sample was extracted using TRIzol reagent (Tiangen, Beijing). The concentration and integrity of RNA were assessed using Qubit 3.0 Fluorometer and Agilent 2100 bioanalyzer. Each group of queen and worker larvae had 3 biological replicates.

Total mRNA of each sample was enriched using oligo (dT) magnetic beads. First-strand cDNA was synthesized by random hexamers, and the second-strand cDNA was synthesized in DNA polymerase I system using dNTPs and RNaseH. Afterwards, cDNA was subjected to end repair and performed with poly(A) tail and ligation sequencing adapter. 200–350 bp cDNA were selected by AMPure XP beads, and then were PCR amplified and purified by AMPure XP beads. Library quality was evaluated using Agilent 2100 bioanalyzer and qRT-PCR. Totally 18 libraries were sequenced by an Illumina NovaSeq 6000 platform. This was performed by Wuhan Benagen Tech Solutions Company Limited also.

TaqMan real-time PCR

The honeybee sample collection method (2d, 4d and 6d queen and worker larvae) was same to that of RNA-Seq, and all samples were from the same colonies of DRS and RNA-Seq. Each larval group had 3 biological replicates, and each biological replicate had 4 technical replicates in the TaqMan real-time PCR. Total RNA of each sample was extracted using TRIzol reagent (Tiangen, Beijing). The cDNA of each sample was synthesized from the total RNA using the Primer-Script RT reagent Kit (TaKaRa) according to the manufacturer’s instructions.

Two isoforms of honeybee Ecr gene (Ecr.t1 and Ecr.t2) were selected for the TaqMan real time PCR verification experiment, and the honeybee GAPDH gene was selected as the house keeping gene. Since the sequences of Ecr.t1 and Ecr.t2 were highly coincident, the TaqMan probes were designed spanning two isoform-specific exons in order to specifically detect these isoforms.The primers were designed that its produces fully covered the probe detecting sequence. The fluorophor and quenching group of probes were 5′-FAM (green, excitation wavelength: 494 nm, emission wavelength: 522) and BHQ-1 respectively. The primers and probes were designed using Primerpremier 5 (version 5.0) and produced by Shanghai GenePharma Co., Ltd (Shanghai, China).

CFX Connect™ real-time PCR (Bio-Rad, USA) was used for this study. Each 20 mL reaction tube contained 10 μL AceQ Universal U+ Probe Master Mix V2, 3 μL cDNA, 0.4 μL Primer Set (10uM), 0.4 μL probe (10uM) and 6.2 μL RNase Free ddH20. The cycling conditions were: 3 min incubation period at 95°C followed by 40 cycles of PCR, each cycle consisting of 12 s at 94°C, 30 s at 59°C and 30 s at 72°C. A cycle threshold (Ct) was calculated by determining the point at which the fluorescence exceeded a threshold limit.

Quantification and statistical analysis

For the data analysis of TaqMan real-time PCR, relative expression of these two isoforms (Ecr.t1 and Ecr.t2) were calculated using 2−ΔΔCt comparative Ct method and were transformed by taking their square root to be normally distributed. Data were analyzed by one-way ANOVA test in Statview package (version 5.0.0.0, SAS institute inc). The p value < 0.05 was considered as significantly different. All groups (2Q, 2W, 4Q, 4W, 6Q and 6W) had three biological replicates from three different honeybee colonies, and each biological replicate had four technical replicates.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Honeybee (Apis mellifera) queen and worker larvae	This paper	N/A
Critical commercial assays
Total RNA Kit I (For DRS)	Omega Bio-tek (USA)	Cat#:R6834
NEBNext Poly(A) mRNA Magnetic Isolation Module	NEB (USA)	Cat#:E7490S
Oxford Nanopore DRS protocol	Oxford Nanopore Technologies (UK)	Cat#:SQK-RNA002
Qubit RNA HS Assay Kit	Thermo Fisher Scientific (USA)	Cat#: Q32852
Reverse Transcriptase Adapter	Oxford Nanopore Technologies (UK)	Cat#:SQK-RNA002
dNTP solution	New England Biolabs (USA)	Cat#:NEB N0447
Quick Ligation Module	New England Biolabs (USA)	Cat#: E6056
Super-Script III Reverse Transcriptase	Thermo Fisher Scientific (USA)	Cat#:18080044
RNA Adapter Mix	Oxford Nanopore Technologies (UK)	Cat#:SQK-RNA002
TRIzol reagent (For RNA-Seq)	Tiangen (China)	Cat#: DP424
oligo (dT) magnetic beads	Vazyme (China)	Cat#: N401-2
RNA Reverse Transcriptase	Thermo Fisher Scientific (USA)	Cat#: 18080044
DNA polymerase I system	New England Biolabs (USA)	Cat#: M0209S
AMPure XP beads	Beckmancoulter (USA)	Cat#A63882
Taq 2X Master Mix	New England Biolabs (USA)	Cat#M0287
Primer-Script RT reagent Kit	TaKaRa (Japan)	Cat#RR037A
AceQ universal U+ probe master mix V2	Vazyme (China)	Cat#Q513-01
50×TAE buffer solution	Solarbio (China)	Cat#T1060
6×DNA Loading Buffer	TransGen (China)	Cat#GH101-01
50 bp DNA Ladder	Tiangen (China)	Cat#MD108
Gsafe Red plus nucleic acid dye	Glpbio (USA)	Cat#GK20002
Deposited data
Direct RNA sequencing data	NCBI database	BioProject: PRJNA748829
RNA-Seq data	NCBI database	BioProject: PRJNA748829
Experimental models: Organisms/strains
Honey bees (Apis mellifera)	Honeybee Research institute, Jiangxi Agricultural university (China)	N/A
Software and algorithms
Guppy	Fesenko et al. (2020)	version 3.2.6
Nanofilt	Coster et al. (2018)	version 2.7.1
Fclmr2	Wang et al. (2018)	version 0.1.2
Minimap2	Li (2018)	version 2.17-r941
Samtools	Li et al. (2009)	version 1.10
Flair collapse	Tang et al. (2020)	version 1.5.0
Stringtie	Pertea et al. (2015)	version 2.1.4
Gffcompare	Pertea and Pertea (2020)	version 0.12.1
Salmon	Patro et al. (2017)	version 1.4.0
DESeq2	Love et al. (2014)	version 1.26.0
KOBAS	Xie et al. (2011)	version 2.0
SUPPA2	Trincado et al. (2018)	version 2.3
Nanopolish	Coster et al. (2018)	version 0.13.2
R package	R Core Team	version 4.0.3
Statview	SAS institute Inc	version 5.0.0.0
PrimerPremier 5	Premier (Canada)	version 5.0
Others
Queen cell frame and excluder	Pan et al., (2013)	N/A
Worker cell frame and excluder	Pan et al., (2013)	N/A